ADAPTIVE NUMERICAL TREATMENT OF 
ELLIPTIC SYSTEMS ON MANIFOLDS 

MICHAEL HOLST 

ABSTRACT. Adaptive multilevel finite element methods are developed and analyzed for 
certain elliptic systems arising in geometric analysis and general relativity. This class of 
nonlinear elliptic systems of tensor equations on manifolds is first reviewed, and then 
adaptive multilevel finite element methods for approximating solutions to this class of 
problems are considered in some detail. Two a posteriori error indicators are derived, 
based on local residuals and on global linearized adjoint or dual problems. The design 
of Manifold Code (MC) is then discussed; MC is an adaptive multilevel finite element 
software package for 2- and 3-manifolds developed over several years at Caltech and UC 
San Diego. It employs a posteriori error estimation, adaptive simplex subdivision, un- 
structured algebraic multilevel methods, global inexact Newton methods, and numerical 
continuation methods for the numerical solution of nonlinear covariant elliptic systems 
on 2- and 3-manifolds. Some of the more interesting features of MC are described in 
detail, including some new ideas for topology and geometry representation in simplex 
meshes, and an unusual partition of unity-based method for exploiting parallel com- 
puters. A short example is then given which involves the Hamiltonian and momentum 
constraints in the Einstein equations, a representative nonlinear 4-component covariant 
elliptic system on a Riemannian 3-manifold which arises in general relativity. A number 
of operator properties and solvability results recently established are first summarized, 
making possible two quasi-optimal a priori error estimates for Galerkin approximations 
which are then derived. These two results complete the theoretical framework for effec- 
tive use of adaptive multilevel finite element methods. A sample calculation using the 
MC software is then presented. 
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1. Introduction 

In this paper we consider adaptive multilevel finite element methods for certain elliptic 
systems arising in geometric analysis and general relativity. Our interest is in developing 
adaptive approximation techniques for the highly accurate and efficient numerical solu- 
tion of this class of problems. We begin by giving a brief introduction to this class of 
nonlinear elliptic systems of tensor equations on manifolds, and then discuss adaptive 
multilevel finite element methods for approximating solutions to this class of problems. 
We derive two a posteriori error indicators, the first of which is local residual-based, 
whereas the second is based on a global linearized adjoint or dual problem. 

The design of a computer program called Manifold Code (MC) is then described, 
which is an adaptive multilevel finite element software package for partial differential 
equations (PDEs) on 2- and 3-manifolds developed over several years at Caltech and UC 
San Diego. MC employs a posteriori error estimation, adaptive simplex subdivision, 
unstructured algebraic multilevel methods, global inexact Newton methods, and numer- 
ical continuation methods for the accurate and efficient numerical solution of nonlinear 
covariant elliptic systems on 2- and 3-manifolds. We describe some of the more inter- 
esting features of MC in detail, including some new ideas for topology and geometry 
representation in simplex meshes. We also describe an unusual partition of unity-based 
method in MC for using parallel computers in an adaptive setting, based on joint work 
with R. Bank [9]. Global L 2 - and iJ 1 -error estimates are derived for solutions produced 
by MC's parallel algorithm by using Babuska and Melenk's Partition of Unity Method 
(PUM) error analysis framework [5] and by exploiting the recent results of Xu and Zhou 
on local error estimation [97]. 

We finish with an example involving the Hamiltonian and momentum constraints in 
the Einstein equations, a representative nonlinear 4-component covariant elliptic system 
on a Riemannian 3-manifold which arises in general relativity. We first summarize a 
number of operator properties and solvability results which were established recently 
in [56]. We then derive two quasi-optimal a priori error estimates for Galerkin approx- 
imations of the constrains, completing the theoretical framework for effective use of 
adaptive multilevel finite element methods. We then present a sample calculation using 
the MC software for this application. More detailed examples involving the use of MC 
for the Einstein constraints may be found in [55, 26]. Applications of MC to problems 
in other areas such as biology and elasticity can be found in [54, 8, 9]. 

2. Adaptive Multilevel Finite Element Methods for Nonlinear 

Elliptic Equations 

In this section we will first give an overview of nonlinear elliptic equations on mani- 
folds, followed by a brief description of adaptive multilevel finite element techniques for 
such equations. 

2.1. Nonlinear elliptic equations on manifolds. Let (Ai, g a b) be a connected compact 
Riemannian rf-manifold with boundary (dAi, cr a b), where the boundary metric o ab is 
inherited from g ab . To allow for general boundary conditions, we will view the boundary 
(d — l)-submanifold dAi (which we assume to be oriented) as being formed from two 
disjoint submanifolds doAi and d\Ai, i.e., 

d M U diM = dAi, d Mr\d!M = 0. (2.1) 

When convenient in the discussions below, one of the two submanifolds d^Ai or d\Ai 
may be allowed to shrink to zero measure, leaving the other to cover dAi . Moreover, 
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in what follows it will usually be necessary to make smoothness assumptions about the 
boundary submanifold dM, such as Lipschitz continuity (for a precise definition see [1]). 
We will employ the abstract index notation (cf. [92]) and summation convention for 
tensor expressions below, with indices running from 1 to d unless otherwise noted. The 
summation convention is that all repeated symbols in products imply a sum over that 
index. Partial differentiation on a non-flat manifold must be covariant, meaning that 
application of gradient and divergence operators require the use of a connection due to 
the curvilinear nature of the coordinate system used to describe the domain manifold. 
Christoffel symbols formed with respect to the given metric g a b (at times denoted %&) 
provide default connection coefficients. 

Covariant partial differentiation of a tensor t a ^,,°£ using the connection provided by the 

metric g a b will be denoted as t a ^"°^ or as D c t a b 1 ^, c b v . Denoting the outward unit normal 
to dM as rib, recall the Divergence Theorem for a vector field w b on M (cf. [65]): 

/ w b . b dx = w b n b ds, (2.2) 
Jm ' JdM 

where dx denotes the measure on M generated by the volume element of g a b\ 

dx = \/ det g a b dx 1 ■ ■ ■ dx d , (2.3) 

and where ds denotes the boundary measure on dM generated by the boundary volume 
element of a ab . Making the choice w b = u ai _ ak v ai - akb in (2.2) and forming the diver- 
gence w b . b by applying the product rule leads to a useful integration-by-parts formula for 
certain contractions of tensors: 

/ u ai ...a k vT akb dx = [ u ai ... ak v"* b n h ds (2.4) 

JM JdM 

- [ V a ^ b U ai ... ak ,b dx. 
JM 

When k = this reduces to the familiar case where u and v are scalars. 

2.1.1. Coupled elliptic systems and augumented systems. Consider now a general second- 
order elliptic system of tensor equations in strong divergence form over M: 

-A ia (x b ,u j ,u k . c ,X)., a + B i (x b ,u j ,u k . c ,X) = OinM, (2.5) 

A ia {x b ,u j ,u k . c ,X)n a + C l {x b ,u j ,u k c ,X) = OondiM (2.6) 

u\x b ) = E\x b ,\) ond M, (2.7) 

where 

AgM™, l<a,6,c<rf, l<i,j,k<n, 

A : M x R n x R nd x R m h-> R nd , B : M x R n x R nd x R m h-> R'\ 

C : d x M x R n x R nd xR m 4K", E : d M xi m 4 R n . 

The divergence-form system (2.5)-(2.7), together with the boundary conditions, can be 
viewed as an operator equation of the form 

G(u, A) = 0, G : Bi x R m h-x B* 2 , (2.8) 

for some Banach spaces B\ and B 2 , where B\ denotes the dual space of B 2 . Analysis 
and numerical techniques often require the Gateaux-linearization operator D u G(u) 6 
C(B\, B 2 ). If D u G(u , X ) is a linear homeomorphism from B\ to B 2 , then the Implicit 
Function Theorem guarantees that there is a neighborhood of (Ao, uq) E R m x B\ con- 
taining regular solutions to (2.8). If D u G(uq, Aq) is singular, so that the Implicit Function 
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Theorem does not apply, then the standard approach is to expand the solution spaces in 
such a way that the expanded problem has a regular solution. This is called the aug- 
mented or bordered system approach to handling folds and bifurcations, and is the basis 
for sophisticated numerical path-following algorithms [62, 63]. In the case of a simple 
limit point (or fold), where D u G{uq) is a Fredholm operator of B\ into B* 2 with index m, 
then the augmented system approach involves simply adding a set of m linear constraints 
to the original system, producing: 



F(u,X,s) 



G{u,\) 
N(u,X,s) 



0, F : Bi x 



M- B* 2 X 



(2.9) 



n i— >• M. m is chosen correctly, then 
for the whole system, which can 



(2.10) 



If the augmentation function N(u, A, s) : B\ x W m x 
the linearization operator DF G C{B% x IR m , x 1 
be written as 

DF(„ X ^ - [ D u G(u,X) D x G(u,X) 
1 ' |_ D«iV(u, A, s) D x N(u, A, s) 

becomes a homeomorphism again. 

Our interest here is primarily in coupled systems of one or more scalar field equations 
and one or more d- vector field equations, possibly augmented as in (2.9)-(2.10). The 
unknown n- vector u % then in general consists of n s scalars and n v d- vectors, so that 
n = n s + n v ■ d. To allow the n-component system (2.5)-(2.7) to be treated notationally 
as if it were a single n-vector equation, it will be convenient to introduce the following 
notation for the unknown vector u % and for the metric of the product space of scalar and 
vector components of u l : 





"«2? 
























9 { 2 ] . 







(2.11) 



If u? k) is a d- vector we take = g a b, if ufa is a scalar we take g^ b ' = 1 



(*.:) 



(*) 



2.1.2. Weak formulations. The weak form of (2.5)-(2.7) is obtained by taking the L 2 - 
based duality pairing between a vector t> J (vanishing on d^M.) lying in a product space 
of scalars and tensors, and the residual of the tensor system (2.5), yielding: 



/ (B* - A ia . a ) v> dx = 0. 
Jm 



(2.12) 



Due to the definition of Gij in (2.11), this is simply a sum of integrals of scalars, each 
of which is a contraction of the type appearing on the left side in (2.4). Using then (2.4) 
and (2.6) together in (2.12), and recalling that v l = on d A4 satisfying (2.1), yields 



(2.13) 



GijA M v 3 . a dx+ / GijB l v 3 dx+ / GijC'v 3 ds = 0. 

M JM JdxM 

Equation (2.13) leads to a covariant weak formulation of the problem: 
Find u e u + B 1 s.t. (F(u), v) = 0, V v G B 2 , 



(2.14) 



for suitable Banach spaces of functions B\ and £>2, where the nonlinear weak form 
(F(-), •) can be written as: 

(F(u), v)= [ Gij(A ia v j . a + BV) dx + [ G %j C i v 3 ds. (2.15) 
Jm JdiM 
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The notation (w, v) will represent the duality pairing of a function v in a Banach space 
B with a bounded linear functional (or form) w in the dual space B*. Depending on the 
particular function spaces involved, the pairing may be thought of as coinciding with 
the L 2 -inner-product through the Riesz Representation Theorem [98]. The affine shift 
tensor u in (2.14) represents the essential or Dirichlet part of the boundary condition if 
there is one; the existence of u such that E = u\d Q M m the sense of the Trace operator is 
guaranteed by the Trace Theorem for Sobolev spaces on manifolds with boundary [93], 
as long as E % in (2.7) and d Ai are smooth enough. If normalization is required as 
in (2.9), then the weak formulation also reflects the normalization. 

2.1.3. Sobolev spaces of tensors. The Banach spaces which arise naturally as solution 
spaces for the class of nonlinear elliptic systems in (2.14) are product spaces of the 
Sobolev spaces Wq This is due to the fact that under suitable growth conditions 

on the nonlinearities in F, it can be shown (essentially by applying the Holder inequality) 
that there exists Pk,Qk, rk satisfying 1 r& < oo such that the choice 

B 1 = W$ x • ■ ■ x W^ D n % B 2 = Wl$ x • ■ ■ x (2.16) 

— + — = 1, r k > mm{p k ,q k }, k = 1, ... ,n e , (2.17) 
Pk Qk 

ensures (F(u),v) in (2.15) remains finite for all arguments [45]. 

The Sobolev spaces are also fundamental to the theory of the finite element method, 
which is based essentially on subspace projection and best approximation. The Sobolev 
spaces W k ' p (M) of tensors on manifolds, and the various subspaces such as W %(M) 
which we will need to make use of later in the paper, can be defined as follows (cf. [51, 
4, 52] for more complete discussions). For a type (r, s)-tensor T a h 1( ^'''^ r = T ! j, where / 
and J are (tensor) multi-indices satisfying |/| = r, \ J\ = s, define 

\T I J \ = (T I J T L M g IL g JM ) 1/2 . (2.18) 
Here, gu and g IJ are generated from the Riemannian d-metric g ab on M. as follows: 

gu = g ab g cd ■ ■ ■ g m , g IJ = g ab g cd ■ ■ ■ g m , (2.19) 

where |/| = \J\ = m, producing m terms in each product. Expression (2.18) is just 
an extension of the Euclidean / 2 -norm for vectors in IR d . For example, in the case of a 
3-manifold, taking |/| = 1, | J\ = 0, g ab = 5 ab , gives: 

\T*j\ = \T a \ = (T a T b g ab ) 1/2 = (T a T% b ) l/2 = \\T a \\ m3) . 

Covariant (distributional) differentiation of order m = \K\ (for some tensor multi- 
index K) using a connection generated by g ab , or generated by possibly a different metric, 
is denoted as any of: 

D m T I J = DkT^ = T 1 ^, (2.20) 

where m should not be confused with a tensor index. Employing the measure dx on M. 
defined in (2.3), the L p -norm of a tensor on M. is defined as: 

WT'jWlhm) = (J iT'jf dx^j ^ , (2.21) 

and the resulting L p -spaces for 1 < p < oo are defined as: 

L p (M) = { T'j | IITSIU.cm) < oo } . (2.22) 

When discussing the properties of //-functions over a manifold Ai we will use the no- 
tation a.e., meaning that the property is understood to hold "almost everywhere" in the 
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sense of Lebesgue measure. We will at times need to make use of the (extended) Holder 
and Minkowski inequalities for tensors in L p -spaces: 

\\U I jW J I \\ L r(M) < ||^j||i>(M)||Wj|U«(A4), (2.23) 

\\U\j + Vj\\ L p( M ) < \\U J jI|lp(.m) + ||Vj||zp(.M), (2.24) 

which hold when U*j, V T j G LP(M), W T j G L q (M), 1/p + l/q = 1/r, 1 < p, q, r < 
oo. The Holder inequality also extends to the case p = 1, q = oo, r = 1, where 
\\ uI j\\l°°(m) = ess sup^g^ \U 7 j(x)|. 
The Sobolev semi-norm of a tensor is defined through (2.21) as: 

\T j\w™, P (m) = ^ J;k\\lp(m)' (2-25) 

|Jf|=m 

and the Sobolev norm is subsequently defined using (2.25) as: 

( i p \ 1/P 

\0<m<fc / 

The resulting Sobolev spaces of tensors are then defined using (2.26) as: 

W k *(M) = {T J j\ WT^Ww^m) < oo } , (2.27) 

Wq' p {M) = { Completion of C^{M) w.r.t. || • \\ W k, P{M) } , (2.28) 

where C^(Ai) is the space of C°° -tensors with compact support in Ai. The space 
W k ' p (M) in (2.28) is a special case of W^{M), which can be characterized as: 

W k ; p {M) = {T j G W k ' p | tr T 7 J; ^ = on d M, \K\ < k - 1} . (2.29) 

Note that if the metric used to define covariant differentiation in (2.20) is taken to be 
different from the metric g a b used in (2.19), it can still be shown that the norms gener- 
ated by (2.26) are equivalent, so that the resulting Sobolev spaces have exactly the same 
topologies [51]. 

The Hilbert space special case of p = 2 is given a simplified notation: 

H k {M) = W k ' 2 {M), (2.30) 

with the same convention used for the various subspaces of H k (Ai) such as Hq(M) and 
Hq d (M). The norm on H k (M) defined above is then actually induced by an inner- 
product as follows: = (T J j, T^Yjjl^y where 

(T / J ,^ J ) L2(A ,)= / T J jS L M g IL g JM dx, (2.31) 

J M 

and where 

i T j> $ I j)Hk{M) = ^J-^k^^j-kIv^m)- (2.32) 

0<\K\<k 

Finally, note that Sobolev trace spaces of tensors living on boundary submanifolds as 
needed for discussing boundary-value problems can be defined under some smoothness 
assumptions on the boundary, and spaces based on fractional-order differentiation (take 
k G E in the discussion above) can be defined in several different ways (cf. [1, 4]). 
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<k{x,y) = l-x-y 
<h,(x,y) = x 
Mx,y) = y 



1 — x — y — z 

y 

x 

z 

Figure 1 . Canonical linear references bases. 

2.2. Adaptive multilevel finite element methods for nonlinear elliptic systems. A 

Petrov-Galerkin approximation of the solution to (2.14) is the solution to the following 
subspace problem: 

Find u h e u h + U h C E x s.t. (F(u h ), v) = 0, V v G V h C B 2 , (2.33) 

for some chosen subspaces Uh and Vh, where dxm(Uh) = dim(V/ l ) — n. A Galerkin 
approximation refers to the case that Uh = Vh- A finite element method is a Petrov- 
Galerkin or Galerkin method in which the subspaces Uh and Vh are chosen to have the 
extremely simple form of continuous piecewise polynomials with local support, defined 
over a disjoint covering of the domain manifold M. by elements. A global C7°-basis 
on the manifold may be defined element-wise from local basis functions defined on a 
reference simplex by use of the chart structure provided with the manifold. For example, 
in the case of continuous piecewise linear polynomials on 2-simplices (triangles) or 3- 
simplices (tetrahedra), the reference element is equipped with the usual basis as shown 
in Figure 2.2. The chart structure provides mappings between the elements contained in 
each coordinate patch and the unit simplex. If the manifold domain can be triangulated 
exactly with simplex elements (possibly as a polyhedral approximation to an underlying 
smooth surface), then the coordinate transformations are simply affine transformations. 
In this sense, finite element methods are by their very nature defined in a chart-wise 
manner. Algorithms for smooth (C h ) 2-surface representations using manifolds have 
been considered recently in [48, 47]; some interesting related work appeared in [37, 38]. 

Due to the non- smooth behavior of their derivatives along simplex vertices, edges, 
and faces in the disjoint simplex covering of M., such continuous piecewise polynomial 
bases clearly do not span a subspace of C l (M); however, one can show [35] that in fact: 

V h = span{0 1 , ...,&,} C Wl${M), M C R d , 

so that continuous, piecewise defined, low-order polynomial spaces do in fact form a 
subspace of the solution space to the weak formulation of the class of second order 
elliptic equations of interest. Making then the choice U h = span{0 1; (f> 2 , . . . , </>„}, V h = 
span{-0i, ip2,..., 4> n }, equation (2.33) in the case of a scalar unknown reduces to a set 
of n nonlinear algebraic relations (implicitly defined) for the n coefficients {ctj} in the 
expansion 

n 

u h = Uh + ^2 a i ( t ) ^ (2.34) 

3=1 

with suitable modification for a vector unknown. In particular, regardless of the com- 
plexity of the form (F(u),v), as long as we can evaluate it for given u and v, then we 





4> 2 (x,y,z) = 
h(x,y,z) = 
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can evaluate the discrete nonlinear residual of the finite element approximation Uh as: 

n 

Fi = (F(u h + ^2a j (f) j ),ip i ), i = 1, ... ,n. 

3=1 

Since the form (F(u),v) involves an integral in this setting, if we employ quadrature then 
we can simply sample the integrand at quadrature points; this is a standard technique in 
finite element technology. Given the local support nature of the functions (fij and ipi, all 
but a small constant number of terms in the sum X}j=i a j ( t ) j are zero at a particular spatial 
point in the domain, so that the residual Fi is inexpensive to evaluate when quadrature is 
employed. 

The two primary issues in using the approximation method are: 

(1) Functionals E[u — Uh) of the error u — Uh (such as norms), and 

(2) Complexity of solving the n nonlinear algebraic equations. 

The first of these issues represents the core of finite element approximation theory, which 
itself rests on the results of classical approximation theory. Classical references to both 
topics include [35, 41, 39]. The second issue is addressed by the complexity theory of 
direct and iterative solution methods for sparse systems of linear and nonlinear algebraic 
equations, cf. [50, 78]. 

2.2.1. Approximation quality: error estimation and adaptive methods. A priori error 
analysis for the finite element method for addressing the first issue is now a very well- 
understood subject [35, 28]. Much activity has recently been centered around a posteriori 
error estimation and the use of error indicators based on such estimates in conjunction 
with adaptive mesh refinement algorithms [10, 7, 6, 90, 91, 97]. These indicators include 
weak and strong residual-based indicators [7, 6, 90], indicators based on the solution of 
local problems [18, 20], and indicators based on the solution of global (but linearized) 
adjoint or dual problems [43]. The challenge for a numerical method is to be as efficient 
as possible, and a posteriori estimates are a basic tool in deciding which parts of the 
solution require additional attention. While the majority of the work on a posteriori 
estimates and indicators has been for linear problems, nonlinear extensions are possible 
through linearization theorems (cf. [90, 91]). The typical solve-estimate -refine structure 
in simplex-based adaptive finite element codes exploiting these a posteriori indicators is 
illustrated in Algorithm 2.2.1. 

Algorithm: (Adaptive multilevel finite element approximation) 

• While (E(u— Uh) is "large") do: 

(1) Find u h E Uh + Uh C. E\ such that (F(uh), v) — 0, V v E Vh C Bi. 

(2) Estimate S(u — Uh) over each element. 

(3) Initialize two temporary simplex lists as empty: Ql = Q2 = 0. 

(4) Simplices which fail an indicator test using equi- distribution of the chosen 
error functional £(u — Uh) are placed on the "refinement" list Ql. 

(5) Bisect all simplices in Ql (removing them from Ql), and place any noncon- 
forming simplices created on the list Q2. 

(6) Ql is now empty; set Ql = Q2, Q2 = 0. 

(7) IfQlis not empty, goto (5). 

• End While. 

The conformity loop (5)-(7), required to produce a globally "conforming" mesh (de- 
scribed below) at the end of a refinement step, is guaranteed to terminate in a finite num- 
ber of steps (cf. [79, 80]), so that the refinements remain local. Element shape is crucial 
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Subdividing 2-simplices (triangles) 



Subdividing 3-simplices (tetrahedra) 




Figure 2. Refinement of 2- and 3-simplices using 4-section, 8-section, 
and bisection. 

for approximation quality; the bisection procedure in step (5) is guaranteed to produce 
nondegenerate families if the longest edge is bisected in two dimensions [81, 88], and if 
marking or homogeneity methods are used in three dimensions [3, 73, 22, 21, 67, 71]. 
Whether longest edge bisection is nondegenerate in three dimensions apparently remains 
an open question. Figure 2 shows a single subdivision of a 2-simplex or a 3-simplex 
using either 4-section (left-most figure), 8-section (fourth figure from the left), or bi- 
section (third figure from the left, and the right-most figure). The paired triangle in the 
2-simplex case of Figure 2 illustrates the nature of conformity and its violation during 
refinement. A globally conforming simplex mesh is defined as a collection of simplices 
which meet only at vertices and faces; for example, removing the dotted bisection in the 
third group from the left in Figure 2 produces a non-conforming mesh. Non-conforming 
simplex meshes create several theoretical as well as practical implementation difficul- 
ties; while the queue-swapping presented in Algorithm 2.2.1 above is a feature unique 
to MC (see Section 3), an equivalent approach is taken in PLTMG [10] and similar 
packages [73, 25, 27, 24]. 



2.2.2. Computational complexity: solving linear and nonlinear systems. Addressing the 
complexity of Algorithm 2.2.1, Newton-like methods as illustrated in Algorithm 2.2.2 are 
often the most effective. 

Algorithm: (Damped-inexact-Newton) 

• Let an initial approximation u be given. 

• While (\(F(u),v) \ > efor any v) do: 

(1) Find w such that (DF(u)w, v) = —(F(u),v) + r, V v. 

(2) Set u = u + \w. 

• End While. 

The bilinear form (DF(u)w, v) which appears in Algorithm 2.2.2 is simply the (Gateaux) 
linearization of the nonlinear form (F(u),v), defined formally as: 



(DF(u)w,v) 



d 
de 



(F(u + ew), v) 



e=0 



This form is easily computed from most nonlinear forms (F(u),v) which arise from sec- 
ond order nonlinear elliptic problems, although the calculation can be tedious in some 
cases (the example we consider later in the paper is in this category). The possibly 
nonzero "residual" term r is to allow for inexactness in the linearization solve for effi- 
ciency, which is quite effective in many cases (cf. [17, 40, 42]). The parameter A brings 
robustness to the algorithm [42, 15, 16]. If folds or bifurcations are present, then the 
iteration is modified to incorporate path-following [62, 14]. 
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As was the case for the Petrov-Galerkin discretized nonlinear residual (F(-), ■), the 
matrix representing the bilinear form in the Newton iteration is easily assembled, regard- 
less of the complexity of the bilinear form (DF(-)-, •). In particular, the matrix equation 

for w = YTj=i Pj4>j nas tne form: 

AU — F, Ui = A, 

where 

n n 

Aij = (DF(u h + 1 ^2a k (f) k )(j) j ,ipi), Fi = (F(u h + ^a j (j) j ),ipi). 

k=l j=l 

As long as the integral-based forms (F(-), ■) and (DF(-)-, ■) can be evaluated at indi- 
vidual points in the domain, then quadrature can be used to build the Newton equations, 
regardless of the complexity of the forms. This is one of the most powerful features of 
the finite element method. It should be noted that there is a subtle difference between the 
approach outlined here (typical for a nonlinear finite element approximation) and that 
usually taken when applying a Newton-iteration to a nonlinear finite difference approxi- 
mation. In particular, in the finite difference setting the discrete equations are linearized 
explicitly by computing the jacobian of the system of nonlinear algebraic equations. In 
the finite element setting, the commutativity of linearization and discretization is ex- 
ploited; the Newton iteration is actually performed in function space, with discretization 
occurring "at the last moment" in Algorithm 2.2.2 above. 

It can be shown that the Newton iteration above is dominated by the computational 
complexity of solving the n linear algebraic equations in each iteration (cf. [17, 49]). 
Multilevel methods are the only known provably optimal or nearly optimal methods for 
solving these types of linear algebraic equations resulting from discretizations of a large 
class of general linear elliptic problems [49, 11, 94]. Unfortunately, the need to accu- 
rately represent complicated PDE coefficient, domain features, and domain boundaries 
with an adapted mesh requires the use of very fine mesh simply to describe the complex- 
ities of the problem, which often precludes the simple solve-estimate-refine approach in 
Algorithm 2.2.1. In Section 3.4 we describe the algebraic multilevel approach we take in 
the MC implementation to adress this, similar to that taken in [31, 32, 83, 89]. 

2.3. Residual-based a posteriori error indicators. There are several approaches to 
adaptive error control, although the approaches based on a posteriori error estimation 
are usually the most effective and most general. While most existing work on a posteri- 
ori estimates has been for linear problems, extensions to the nonlinear case can be made 
through linearization. For example, consider the nonlinear problem in (2.9), which we 
will write as follows (ignoring the parameters for simplicity): 

F(u) = 0, F E C\Bi, B* 2 ), Bi, B 2 Banach spaces, (2.35) 

and a discretization: 

F h {u h ) = 0, F h G C°(U h , V£), U h cB h V h C B 2 . (2.36) 

The nonlinear residual F(uh) can be used to estimate the error \\u — Uh\\Bi> through the 
use of a linearization theorem [68, 90]. An example of such a theorem due to Verfiirth is 
the following. 

Theorem 2.1. [90] Let u e X be a regular solution of F(u) = 0, so that the Gateaux 
derivative DF(u) is a linear homeomorphism of B\ onto B\. Assume DF is Lipschitz 
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continuous at u, so that there exists R such that 



7 



sup 

u h GB(u,R ) 



DF{u)-DF{u h )\\ c{ B u Bi) 



< oo. 



U - Uh\\B 



Let R = min{i? ,7 1 \\DF(u) l \\c{Bi,Br)^l 1 \\DF(u)\\c(b u b*)}- Then for all u h e 
B(u,R), 



The effect of linearization is swept under the rug somewhat by the choice of R suffi- 
ciently small, where R is the radius of an open ball in B\ about it, denoted as B(u, R) 
in the theorem above. One then ignores the factors in (2.37) involving the linearization 
DF(u) and its inverse, and focuses on two-sided estimates for the nonlinear residual 
ll-P 1 ^)!!^* appearing on each side of (2.37). Since one typically constructs highly re- 
fined meshes where needed, such local linearized estimates are thought to reasonable, 
although much evidence to the contrary has been assembled by the dual-problem error 
indicator community (see the discussion later in this section). Note that ||.F(w/i)||b* can 
be estimated in different ways, including 

(1) Approximation by H-F/i^/JHs* (residual estimates) [90, 91, 68, 69], 

(2) Solution of local Neumann (or Dirichlet) problems [18, 20]. 

The approaches can be shown to be essentially equivalent (up to constants; cf. [90, 23]). 
For reasons of efficiency, estimation by strong residuals is often used rather than the 
solution of local problems in the case of elliptic systems and/or in the setting of three- 
dimensional problems. In particular, one employs the linearization theorem above, to- 
gether with some derived (and computable) upper and lower bounds on the nonlinear 
residual ||-F , («/ v )||b| given by the following pair of inequalities: 



While it is clear that the upper bound C4 is the key to bounding the error, the lower bound 
C3 can also be quite useful; it can help to ensure that the adaptive procedure doesn't do 
too much work by over-refining an area where it is unnecessary. The effectiveness of 
an adaptive finite element code can hinge on the implementation details of the estimator, 
and implementing it efficiently can be quite an art form (cf. [90, 18, 20]). 

We now consider the first two of these approaches in more detail. First, we derive 
a strong residual-based a posteriori error indicator for general Petrov-Galerkin approx- 
imations (2.33) to the solutions of general nonlinear elliptic systems of tensors of the 
form (2.5)-(2.7). The analysis involves primarily the weak formulation (2. 14)— (2. 15). 
Our derivation follows closely that of Verfiirth [90, 91] in the flat, Cartesian case. In 
the next section, we will consider the second alternative, namely an indicator based on a 
duality approach. 

It should be noted that while the discussions throughout the paper are generally valid 
for domains which are connected compact Riemannian manifolds with Lipschitz con- 
tinuous boundaries, several of the results we will need to employ here have only been 
shown to hold in the case of bounded 2- and 3-manifolds with smooth boundaries, with 
an atlas consisting of only one chart (i.e., bounded open subsets of IR 2 and IR 3 ). Examples 
are convex polyhedra in W 1 , which automatically satisfy the Lipschitz continuity assump- 
tion. The extensions of some of these results from open sets to Riemannian manifolds are 
not immediate; a number of subtle issues arise when manifold domains are considered in 




'2 ■ 



(2.37) 



□ 



C 3 < \\F{u h )\\ B » 2 <C 4 . 
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conjunction with Sobolev spaces, the subject of two recent monographs [4, 52]. Some of 
the difficulties encountered impact approximation theory on manifolds [1,4, 82, 86, 52]. 
However, if a sufficient amount of the function space framework is in place, and if the 
underlying manifold admits a partition of unity (e.g., if it is paracompact), then it should 
be possible to extend finite element approximation theory by operating chartwise, and 
then globalizing the local results using the partition of unity. We will assume here that 
such extensions are possible; understanding the manifold case of multiple charts and 
other complications is work in progress [53]. 

The starting point for our residual-based error indicator is the linearization inequal- 
ity (2.37). In our setting of the weak formulation (2.14)— (2.15), we make the appropriate 
choice (2. 16)— (2. 17), where we restrict our discussion here to a single elliptic system for 
a scalar or a d-vector (i.e., the product space has dimension n e = 1), which includes the 
examples presented later in the paper. The linearization inequality then involves standard 
Sobolev norms: 

Ci\\F(uh)\\w- l >9(M) < \\ u - Uh\\w^(M) < C 2 \\F(uh)\\w- 1 ><i(M), (2.38) 

for 1/p + 1/q — 1, r > mm{p, q}, where W~ l ' q (M) = (W l,q (M))* denotes the dual 
space of bounded linear functionals on W l,q (M). The norm of the nonlinear residual 
F(-) in the dual space of bounded linear functionals on W 1,q (M) is defined in the usual 
way: 

\\F{u)\\ w -,, q{M) = sup 1 } F(U) ' V)1 . (2.39) 

The numerator is the nonlinear weak form (F(u),v) appearing in (2.15). (We will con- 
sider only the case of no parameters; see [90] for the case of parameters.) In order to 
derive a bound on the weak form in the numerator we must first introduce quite a bit of 
notation that we have managed to avoid until now. 

To begin, we assume that the d-manifold Ai has been exactly triangulated with a 
set S of shape-regular d-simplices (the finite dimension d is arbitrary throughout this 
discussion). A family of simplices will be referred to here as shape-regular if for all 
simplices in the family the ratio of the diameter of the circumscribing sphere to that 
of the inscribing sphere is bounded by an absolute fixed constant, independent of the 
numbers and sizes of the simplices that may be generated through refinements. (For a 
more careful definition of shape-regularity and related concepts, see [35].) It will be 
convenient to introduce the following notation: 

S = Set of shape-regular simplices triangulating M 

Af(s) = Union of faces in simplex set s lying on dj^Ad 

Z(s) = Union of faces in simplex set s not in Af(s) 

F(s) = Af( s )UX(s) 

°°s = [j{sES\sC]s^0, where s € S } 

Uf = \J{§ eS \ / D s 7^ 0, where / e JF } 

h s = Diameter (inscribing sphere) of the simplex s 

hf = Diameter (inscribing sphere) of the face /. 

When the argument to one of the face set functions M, X, or F is in fact the entire set of 
simplices S, we will leave off the explicit dependence on S without danger of confusion. 
Referring forward briefly to Figure 3 will be convenient. The two darkened triangles in 
the left picture in Figure 3 represents the set Wf for the face / shared by the two triangles. 
The clear triangles in the right picture in Figure 3 represents the set w s for the darkened 
triangle s in the center (the set w s also includes the darkened triangle). 



ADAPTIVE NUMERICAL TREATMENT OF ELLIPTIC SYSTEMS ON MANIFOLDS 13 

Finally, we will also need some notation to represent discontinuous jumps in function 
values across faces interior to the triangulation. To begin, for any face / 6 J\f, let n / 
denote the unit outward normal; for any face / G X, take n/ to be an arbitrary (but fixed) 
choice of one of the two possible face normal orientations. Now, for any v E L 2 (A4) 
such that v £ C°(s) Vs E S, define the jump function: 

[v]f(x) = lim v{x + erif) — lim v(x — en/). 

We now begin the analysis by splitting the volume and surface integrals in (2.15) 
into sums of integrals over the individual elements and faces, and we then employ the 
divergence theorem (2.4) to work backward towards the strong form in each element: 



(F(u),v) = [ Q lJ (A i V. a + BV)dx+ [ QijCV ds 

JM Jd N M 

= E / c -'o(- r ' J ,, + £V) dx + Y [ SijCTv 1 ds 

= E f ~ dX + E f Qi^V ds 

seS J s seS ^ 9s 

+ E / g ^° iv3 ds - 

Using the fact that (2.33) holds for the solution to the discrete problem, we employ the 
jump function and write 

(F(u h ),v) = (F(u h ),v-v h ) (2.40) 

= E f SijiB' - A%)(ui - vi) dx 
ses Js 

+ E / Q v A ia n a {v 3 -vi)ds 

+ E / GijC^-vi) ds 
feAf f 

= E [^(B l -A™)(vi-vi)dx 
ses Js 

+ E / GiA^ndf^-vl) ds 
fex J f 

+ E / + ^X)(^' - v{) ds 

< Edi^-^iu^ii^-^iu^)) 

seS 

+ E(ll [^J/IUn/)il^'-4lU"(/)) 

/ex 

+ E(ll^ + ^XIU,(/)lk i -^IU,(/ ) ), 

feAf 

where we have applied the Holder inequality (2.23) three times with 1/p + 1/q = 1. 
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In order to bound the sums on the right, we will employ a standard tool known as 
a H /1,p -quasi-interpolant I h . An example of such an interpolant is due to Scott and 
Zhang [87], which we refer to as the SZ- interpolant (see also Clement's interpolant 
in [36]). Unlike point-wise polynomial interpolation, which is not well-defined for func- 
tions in W^ P (M) when the embedding W 1>P (M) ^ C°(M) fails, the SZ-interpolant 
Ih can be constructed quite generally for VT 1,p -functions on shape-regular meshes of 
2- and 3-simplices. Moreover, it can be shown to have the following remarkable local 
approximation properties: For all v E W 1,q (M), it holds that 

\\v - Ihv\\Li( s ) < C s h s \\v\\ w i, q{uls ), (2.41) 
\\v - I h v\\ Lq{f) < Cfh 1 f 1/9 \\v\\ W i, q ( )ilf y (2.42) 

For the construction of the SZ-interpolant, and for a proof of the approximation inequal- 
ities in L p -spaces for p ^ 2, see [87]. A simple construction and the proof of the first 
inequality can also be found in the appendix of [57]. 

Employing now the SZ-interpolant by taking Vh = IhV in (2.40), using (2.41)-(2.42), 
and noting that 1 — 1/q = 1/p, we have 



(F(u h ),v) < J2 C M Bi ~ A %\\ 



s<=S 



/ex 



+ CfrfV + Aia na\\LPtf)\\vi\\ w ^ } ) 



< 



\Y, CP s hP s\\ Bi ~ A 1a\\l, { s) (2-43) 



+J2° p f h f\\ i Aia ^] f \\Uf) + E c f h f\\ ci + A * n *\ 



/ex feM 




1/9 



/ex feAf 



where we have used the discrete Holder inequality to obtain the last inequality. 

It is not difficult to show (cf. [90]) that the simplex shape regularity assumption bounds 
the number of possible overlaps of the sets cu s with each other, and also bounds the 
number of possible overlaps of the sets cjf with each other. This makes it possible to 
establish the following two inequalities: 

E IMIw^K) - -^IMIw^CM)' (2.44) 

Ell^ll^K) < D fW v \\w^(My ( 2 - 45 ) 
/e-F 

where D s and Df depend on the shape regularity constants reflecting these overlap 
bounds. Therefore, since X C T and J\f C J 7 , we employ (2.44)-(2.45) in (2.43) which 
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l/p 



gives 

(F(u h ),v) < C 5 \\ V \\ W i, HM) ■ [J2 hP s\\ Bi ~ A 1a\\Us) ( 2 - 46 ) 

/ex feAf J 

where C 5 = ma.x SyJ r{C s , C/} • max^^jDy 9 , D 1 / 9 } depends on the shape regularity of 
the simplices in S. 

We finally now use (2.46) in (2.39) to achieve the upper bound in (2.38): 

\\ u — Uh\\w 1 ' r (M) < C 2 \\F(u h )\\ w -i,q( M ) 

— C-2 SUp -r—r 

Q^v^W 1 'i{M) \\ v \\W 1 -i{M) 

< c 2 c 5 (j2 hP s\\ Bi - A %\\Us) ^ 



ip 

1/ Hlp(/) 

/ex 



l/p 



We will make one final transformation that will turn this into a sum of element-wise 
error indicators that will be easier to work with in an implementation. We only need to 
account for the interior face integrals (which would otherwise be counted twice) when 
we combine the sum over the faces into the sum over the elements. This leave us with 
the following 

Theorem 2.2. Let u e W 1,r (Ai) be a regular solution of (2.5)-(2.7), or equivalently 
of (2.14)-(2.15), where (2.16)-(2.17) holds. Then under the same assumptions as in 
Theorem 2.1, the following a posteriori error estimate holds for a Petrov-Galerkin ap- 
proximation Uh satisfying (2.33): 

\\u-u h \\ w u r{M) <c(y2A , (2.48) 



where 

C = 2 • max{C s , C f } ■ m&x{Dl /q , D 1 / 9 } ■ WDF^Waw-^w^ 
and where the element-wise error indicator r\ s is defined as: 



B l - A%\\% {s) + i £ h f \\ [A ia n a ] f \\l P(f) (2.49) 



2 

fex(s) 



i/p 

a\\ L p(f) 



/eA/"(s) 

Proof. The proof follows from (2.47) and the discussion above. □ □ 
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The element-wise error indicator in (2.49) provides an error bound in the W /1 ' r -norm 
for a general covariant nonlinear elliptic system of the form (2.5)-(2.7), with 1 fp+ l/q — 
1, r > min{p, q}, which may be more appropriate than the r = p = q = 2 case for some 
nonlinear problems. Issues related to this topic are discussed in [90]. Following [90], it 
is possible to use a similar analysis to construct lower bounds, dual to (2.48), of the form 



\SScS / 

Such results are useful for performing unrefinement and in accessing the quality of an 
error-indicator. 

2.4. Duality-based a posteriori error indicators. We now derive an alternative a pos- 
teriori error indicator for general Petrov-Galerkin approximations (2.33) to the solutions 
of general nonlinear elliptic systems of tensors of the form (2.5)-(2.7). The indica- 
tor is based on the solution of a global linearized adjoint (or dual) problem; again the 
analysis involves primarily the weak formulation (2.14)— (2. 15), and follows closely that 
of [23, 43]. This approach can be viewed as simply another way to bound the nonlinear 
residual after employing the (possibly quite crude) one-time linearization in 

Theorem 2.1. However, the approach can be used to avoid the one-time linearization 
step, bringing the stability properties of the differential operator into the error indica- 
tor by updating the linearization as the solution is improved, and by incorporating the 
linearization operator itself into the error indicator. 

As before, we are interested in the solution to the operator equation (2.35) and also in 
error estimates for approximations Uh satisfying (2.36). We begin with the generalized 
Taylor remainder in integral form: 

F{u + h) = F{u) + | jf DF(u + £/i)<f£ j h. (2.50) 

Taking h = Uh — u, the error e = u — Uh can be expressed as follows: 

R = -F(u h ) = -F(u + [u h - u]) = -F(u) - A(u h - u) = - Ae, 
where the linearization operator A is defined from (2.50) as: 



If a linear functional of the error 1(e) = (e, ip) is of interest rather than the error itself, 
where ij) is the Riesz-representer of /(•), then we can exploit the linearization operator A 
in (2.51), and its (unique) adjoint A T , to produce an error indicator: 



for the residual weights <p, where the data for the dual problem ip is the Riesz-representer 
of the functional of interest. Strong norm estimates of the form (2.38) can be established 
using duality (cf. [23]), but the operator information represented by the dual solution cf> 
is then lost (it appears in the constants). If a functional of the error is of interest (e.g., the 
error along a curve or surface in the domain), then a more delicate approach is to instead 
employ the dual solution as part of the indicator: 





(2.51) 



|(e,V}| = \(e,A T 4>)\ = \(Ae,4>)\ = \(R,<f>)\ = \(F(u h ),<f>)\. 
The indicator requires the solution of the linearized dual problem: 

A T (j) = i) 



(2.52) 



|(e, ip) \ = \(F(uh), 4>)\ < error estimate. 



(2.53) 
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The dual solution <\> obtained by solving (2.52) is used locally (element-wise) in (2.53), 
with the dual solution as residual weights (cf. [43]). 

To construct such estimates for general Petrov-Galerkin approximations (2.33) to the 
solutions of general nonlinear elliptic systems of tensors of the form (2.5)-(2.7), we first 
need some simple identities to help identify the form of the linearized dual problem: 

A ia {u k ,u k c ) - A ia {U k ,U k c ) 

1 d 

—A ia (su k + (1 - s)U k } su k c + (1 - s)U k c )ds 
q ds 

[ {D 1 A ia {su k + (1 - s)U k , su k . c + (1 - s)U k c ) 
Jo 

~[su k + (1 - s)U k ] 

+D 2 A ia {su k + (1 - s)U k ,su k . c + (1 - s)U k c ) 
~[«t* c + (l-«)?7*j}cfa 

jjT D 1 A ia (su k + (1 - s)U k , su k c + (1 - s)U k c )ds^ {u k - U k ) 
+ | jT D 2 A ia (su k + (1 - s)U k , su k . c + (1 - s)U k c )ds^ {u k c - U k ) 



A ia b e b + A% c e b . c . 



B l (u k ,u k . c )- B l (U k ,U k c ) 

1 d 

j- s B\su k + (1 - s)U k , su k . c + (1 - s)U k c )ds 



[ {DxB^su* + {l-s)U k , su k . c + (1 - s)U k c ) 
Jo 

~[su k + (1 - s)U k ] 

+D 2 B\su k + (1 - s)U k , su k . c + (1 - s)U%) 
~ s [su) c + (l-s)U%]\ds 

j jf D 1 B i (su k + (1 - s)U k , su k c + (1 - s^jdsj - U k ) 
i ^ 



+ |y D 2 B\su k + (i - s)c/ fc , s< c + (i - s)c/^ c )ds ;> « - c/; 

Similarly, 

&(u k ) - t7 l '(C/ fc ) = |2 DiC*(s«* + (1 - s)f/ fe )ds }> (w fc - £/ fc ) = CV 6 
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Therefore, given our original weak form in (2.15), we have 

(F(u)-F(U),cf>) = [ Gid[A ia {u\u\ c )-[A ia {U\U\ c )W. a 

J M 

+ [B\u\u\ c )-B\U\U\ c )W) dx 

+ / g VJ \c i {u k )-c\u k )w ds 

JdxM 



+ (B\e b + B\ c e b . c W} dx + f G l3 C\e h & ds 

JdiM 

= (Ae,<f>) 
= (e, A T (f)). 

The weak form of the linearized dual problem is then: 

Find G B 1 such that (A T <f), v) = (if), v), \fv G B 2 , (2.54) 
where the adjoint form is 

(A T (fi, v)= f Ga {A% % a v b . c + AW, a v b + B\ c 4Pv b . c (2.55) 
+ B\4> j v b } dx + [ G l f l h &v h ds. 

Jd x M 

The strong form of the linearized dual problem in (2.54)-(2.55) is then: 

Qii {AW* - (^1 C ^';a) ;c + & b </>> - (W). e } = in M, 

g t] {A ia b c ^; a nc+{B\ c n c + C\)^} = OonfrM, 

u\x h ) = Oond M. 

This leads to the following error representation: 

Theorem 2.3. Given a projector Ph : Si 4 Uh onto the finite element subspace Uh C 
B\, the functional error is: 

where 

(K(U), v) = ^ / ^(^ - A*)^' - j^) (2.56) 



/ex 



feAf f 
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Proof. We begin by working backward toward the strong form: 

(e,V> = (F(U),<j>) 

= I t;,;(.l"'o' :( , • B : <>•)<!.,■ ■ [ Q t] C 1 ^ ds 

JM Jd N M 

= Y, f &A Bi - A%W dx + Y I Gij^n^ ds 
ses J s ses ^ 9s 

+ Y [ g v Gi< P ds - 

Using then Galerkin orthogonality and a jump function gives: 

(F{U),<f>) = (F(U),<f>-P h <f>) 

= E h^ ,r - A*£tf - Prf) dx 

ses Js 

+ Y j GijA ia n a (<fj - P h (fP) ds 

seS Jds 

+ Y [ GifiW - Phft) ds 

= E / g ^ Bi - A v^ J - dx 

ses Js 

+ Y I % [A ia n a ] f {^ -P h 4P)ds 
fez J f 

+ Y [ G ^° l + - Ph¥) ds 

= (n{u),<f>). 

□ □ 

The error representation can be used as an error indicator as illustrated in Algo- 
rithm 2.4. 

Algorithm: (Linearized dual error indicator) 

(1) Decide which linear functional s) of the error 1(e) = (e,ip) is of interest. 

(2) Pose and solve the linearized dual problem A T cj) = ipfor the dual weight function 
0. 

(3) Numerically approximate (TZ(U),<p) within each element as an element-wise er- 
ror indicator. 

(4) Elements which fail an indicator test (using equi-distribution) are marked for 
refinement. 

3. Manifold Code (MC): adaptive multilevel finite element methods 

ON MANIFOLDS 

MC(see also [54, 8, 9, 55, 26]) is an adaptive multilevel finite element software pack- 
age, written in ANSI C, which was developed by the author over several years at Caltech 
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and UC San Diego. It is designed to produce highly accurate numerical solutions to non- 
linear covariant elliptic systems of tensor equations on 2- and 3-manifolds in an optimal 
or nearly-optimal way. MC employs a posteriori error estimation, adaptive simplex sub- 
division, unstructured algebraic multilevel methods, global inexact Newton methods, and 
numerical continuation methods for the highly accurate numerical solution of nonlinear 
covariant elliptic systems on (Riemannian) 2- and 3-manifolds. 

3.1. The overall design of MC. MC is an implementation of Algorithm 2.2.1, where 
Algorithm 2.2.2 is employed for solving nonlinear elliptic systems that arise in Step 1 
of Algorithm 2.2.1. The linear Newton equations in each iteration of Algorithm 2.2.2 
are solved with algebraic multilevel methods, and the algorithm is supplemented with 
a continuation technique when necessary. Several of the features of MC are somewhat 
unusual, allowing for the treatment of very general nonlinear elliptic systems of tensor 
equations on domains with the structure of 2- and 3-manifolds. In particular, some of 
these features are: 

• Abstraction of the elliptic system: The elliptic system is defined only through 
a nonlinear weak form over the domain manifold, along with an associated lin- 
earization form, also defined everywhere on the domain manifold (precisely the 
forms (F(u),v) and (DF(u)w,v) in the discussions above). To use the a pos- 
teriori error indicators, a third function F(u) must also be provided (essentially 
the strong form of the problem). 

• Abstraction of the domain manifold: The domain manifold is specified by giving 
a polyhedral representation of the topology, along with an abstract set of coordi- 
nate labels of the user's interpretation, possibly consisting of multiple charts. MC 
works only with the topology of the domain, the connectivity of the polyhedral 
representation. The geometry of the domain manifold is provided only through 
the form definitions, which contain the manifold metric information, and through 
a oneChart ( ) routine that the user provides to resolve chart boundaries. 

• Dimension independence: Exactly the same code paths in MC are taken for both 
two- and three-dimensional problems (as well as for higher-dimensional prob- 
lems). To achieve this dimension independence, MC employs the simplex as its 
fundamental geometrical object for defining finite element bases. 

As a consequence of the abstract weak form approach to defining the problem, the com- 
plete definition of a complex nonlinear tensor system such as large deformation nonlinear 
elasticity requires writing only a few hundred lines of C to define the two weak forms, 
and to define the oneChart ( ) routine. Changing to a different tensor system (e.g. the 
example later in the paper involving the constraints in the Einstein equations) involves 
providing only a different definition of the forms and a different domain description. 

3.2. Topology and geometry representation in MC: The Ringed Vertex. A datas- 
tructure referred to as the ringed-vertex (cf. [53]) is used to represent meshes of d- 
simplices of arbitrary topology. This datastructure is illustrated in Figure 3. The ringed- 
vertex datastructure is similar to the winged-edge, quad-edge, and edge-facet datas- 
tructures commonly used in the computational geometry community for representing 
2-manifolds [72], but it can be used more generally to represent arbitrary d-manifolds, 
d > 2. It maintains a mesh of <i-simplices with near minimal storage, yet for shape- 
regular (non-degenerate) meshes, it provides 0(l)-time access to all information neces- 
sary for refinement, un-refinement, and Petrov-Galerkin discretization of a differential 
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Figure 3. Polyhedral manifold representation. The figure on the left 
shows two overlapping polyhedral (vertex) charts consisting of the two 
rings of simplices around two vertices sharing an edge. The region con- 
sisting of the two darkened triangles around the face / is denoted Uf, 
and represents the overlap of the two vertex charts. Polyhedral mani- 
fold topology is represented by MC using the ringed-vertex (or RIVER) 
datastructure. The datastructure is illustrated for a given simplex s in the 
figure on the right; the topology primitives are vertices and ci-simplices. 
The collection of the simplices which meet the simplex s at its vertices 
(which then includes those simplices that share faces as well) is denoted 
as cu s . (The seto; s includes s itself.) Edges are temporarily created during 
subdivision but are then destroyed (a similar ring datastructure is used to 
represent the edge topology). 

operator. The ringed-vertex datastructure also allows for dimension independent imple- 
mentations of mesh refinement and mesh manipulation, with one implementation (the 
same code path) covering arbitrary dimension d. An interesting feature of this datastruc- 
ture is that the C structures used for vertices, simplices, and edges are all of fixed size, so 
that a fast array-based implementation is possible, as opposed to a less -efficient list-based 
approach commonly taken for finite element implementations on unstructured meshes. A 
detailed description of the ringed-vertex datastructure, along with a complexity analysis 
of various traversal algorithms, can be found in [53]. 

Since MC is based entirely on the rf-simplex, for adaptive refinement it employs sim- 
plex bisection, using one of the simplex bisection strategies outlined earlier. Bisection 
is first used to refine an initial subset of the simplices in the mesh (selected accord- 
ing to some error indicator combined with equi-distribution, discussed below), and then 
a closure algorithm is performed in which bisection is used recursively on any non- 
conforming simplices, until a conforming mesh is obtained. If it is necessary to improve 
element shape, MC attempts to optimize the following simplex shape measure function 
for a given d-simplex s, in an iterative fashion, similar to the approach taken in [19]: 



7/0, eft = 



2 2 ( 1 -a)3^ 


\s 


2 
d 


2~2o<i<j<d 




2 



(3.1) 



The quantity \s\ represents the (possibly negative) volume of the <i-simplex, and | e^- 1 
represents the length of the edge that connects vertex i to vertex j in the simplex. For 
d = 2 this is the shape-measure used in [19] with a slightly different normalization. For 
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d = 3, the measure in (3.1) is the shape-measure developed in [66] again with a slightly 
different normalization. The shape measure above can be shown to be equivalent to the 
sphere ratio shape measure commonly used (cf. [66]). 

3.3. Discretization, adaptivity, and error estimation in MC. Given a nonlinear weak 
form (F(u),v), its linearization bilinear form (DF(u)w, v), a Dirichlet function u, and 
a collection of simplices representing the domain, MC uses a default linear element to 
produce and then solve the implicitly defined nonlinear algebraic equations for the ba- 
sis function coefficients in the expansion (2.34). The user can also provide their own 
element, specifying the number of degrees of freedom to be located on vertices, edges, 
faces, and in the interior of simplices, along with a quadrature rule, and the values of 
the trial (basis) and test functions at the quadrature points on the master element. Dif- 
ferent element types may be used for different components of a coupled elliptic system. 
The availability of a user-defined general element makes it possible to, for example, use 
quadratic elements as would be required in elasticity applications to avoid locking. 

Once the equations are assembled and solved (discussed below), a posteriori error es- 
timates are computed from the discrete solution to drive adaptive mesh refinement. The 
idea of adaptive error control in finite element methods is to estimate the behavior of the 
actual solution to the problem using only a previously computed numerical solution, and 
then use the estimate to build an improved numerical solution by upping the polynomial 
order (p-refinement) or refining the mesh (/i-refinement) where appropriate. Note that 
this approach to adapting the mesh (or polynomial order) to the local solution behavior 
affects not only approximation quality, but also solution complexity: if a target solution 
accuracy can be obtained with fewer mesh points by their judicious placement in the 
domain, the cost of solving the discrete equations is reduced (sometimes dramatically) 
because the number of unknowns is reduced (again, sometimes dramatically). Gener- 
ally speaking, if an elliptic equation has a solution with local singular behavior, such as 
would result from the presence of abrupt changes in the coefficients of the equation, or 
a domain singularity, then adaptive methods tend to give dramatic improvements over 
non-adaptive methods in terms of accuracy achieved for a given complexity price. Two 
examples illustrating bisection-based adaptivity patterns (driven by a completely geo- 
metrical "error" indicator simply for illustration) are shown in Figure 4. 

MC employs the error indicators derived in Section 2.3 adaptive solution of nonlinear 
elliptic systems of the form (2.5)-(2.7). In particular, the indicators are used to adap- 
tively construct Galerkin solutions satisfying (2.33), which approximate weak solutions 
satisfying (2.14). MC can be directed to use either the local residual indicator (2.49) to- 
gether with the principle of error equi-distribution, or the duality -based weighted residual 
indicator (2.56), again together with equi-distribution. 

Of course, we can't perform the integrals in (2.49) or (2.56) exactly in most cases, so 
we employ quadrature in MC. Another option is to project the data onto the finite ele- 
ment spaces involved and then to perform the integrals exactly; this approach is analyzed 
carefully in [90]. Note that r) s is computable by quadrature, since all terms appearing in 
the definition depend only on the (available) computed solution Uh- In particular, each of 
the terms 

A ia ( x \ (u h y, K) fc ;c ) ;a , A ia ( x \ ( Uh y, K) fc ; > a , 
b\x\ (u h y, K) fc ; j, c\ x \ ( Uh y, M fc ;c ), 

depend only on Uh, its first derivatives, and the normal vector n q (which is known from 
simple geometrical calculations). All of the terms but the first are already provided by 
the user as part of the weak form (F(u),v) required to use MC. The only problematic 
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Figure 4. Examples illustrating the 2D and 3D adaptive mesh refine- 
ment algorithms in MC. The right-most figure in each row shows a close- 
up of the area where most of the refinement occured in each example. 



term is the first one; this represents the strong form of the principle part of the equation, 
and must be supplied by the user as a separate piece of information. 

In order to understand this more completely, we will briefly describe some of the prob- 
lem specification details in MC. To use MC to discretize problems of the form (2.14)- 
(2.15), the user is expected to provide the Dirichlet function u by providing the function 
E in (2.7). In addition, the user provides the nonlinear weak form: 

(F(u),v) = [ g ij (A ia v j . a + B i v j )dx + [ QijCVds 

JM Jd N M 

F (u)(v)dx + F 1 (u)(v)ds, (3.2) 

M Jd N M 

by providing the integrand function F t (u)(v) defined as: 

F t (u)(v) = { G^,-irv% ift = o, 

* WW \ Gy&v 3 , if t = 1. { ' 

In order to use the inexact Newton iteration in MC to produce a Petrov-Galerkin approxi- 
mation satisfying (2.33), the user must also provide a corresponding bilinear linearization 
form: 



(DF(u)w,v) = / DF (u)(w,v) dx + / DF^u^w, v) ds, (3.4) 

JM Jd N M 

where the integrand function DF t (u)(w, v) is defined through Gateaux differentiation as 
described in Section 2.1. In order to use the a posteriori error estimator in MC, the user 
must provide an additional vector- valued function SF t (u), defined as: 

B i - A ia a , if t = 0, 
SF t (u) = { O + A&na, if t = 1, (3.5) 
A ia n a , if t = 2. 
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The key point that must be emphasized here is the following: since MC employs quad- 
rature to evaluate the integrals appearing in each of (3.2), (3.4), and (2.49), the user- 
provided functions Ft(u)(v), DFt(u)(w, v), and SF t (u) only need to be evaluated at a 
single point x p E M at a time. In other words, the user can simply evaluate the ex- 
pressions in (3.3) and in (3.5) as if they were point vectors and point tensors, rather than 
vector and tensor fields. This is one of the most powerful features of MC, and of nonlin- 
ear finite element software in general. It implies that the user-defined functions F t (u)(v), 
DF t (u)(w, v), and SF t (u) can usually be implemented to appear in software exactly as 
they do on paper. 

The remaining quantities appearing in the estimator (2.49), namely the normal vector 
n q and the mesh parameters h s and hf, are completely geometrical and can be computed 
from the local simplex geometry information. The indicator is very inexpensive when 
compared to the typical cost of producing the discrete solution Uh itself; the number 
of function evaluations and arithmetic operations (for performing quadrature) is always 
linear in the total number of simplices. Moreover, the indicator is completely local; it can 
be computed chart-wise when multiple coordinate systems are employed. These ideas 
are explored more fully in [53]. 

3.4. Solution of linear and nonlinear systems in MC. When a system of nonlin- 
ear finite element equations must be solved in MC, the global inexact-Newton Algo- 
rithm 2.2.2 is employed, where the linearization systems are solved by linear multilevel 
methods. When necessary, the Newton procedure in Algorithm 2.2.2 is supplemented 
with a user-defined normalization equation for performing an augmented system contin- 
uation algorithm. The linear systems arising as the Newton equations in each iteration 
of Algorithm 2.2.2 are solved using a completely algebraic multilevel algorithm. Either 
refinement-generated prolongation matrices P k , or user-defined prolongation matrices 
P k in a standard YSMP-row-wise sparse matrix format, are used to define the multilevel 
hierarchy algebraically. In particular, once the single "fine" mesh is used to produce the 
discrete nonlinear problem F{u) = along with its linearization Au = f for use in the 
Newton iteration in Algorithm 2.2.2, a J-level hierarchy of linear problems is produced 
algebraically using the following recursion: 

A k+1 = PlA k P k) k = l,...,J — l, A X = A. 

As a result, the underlying multilevel algorithm is provably convergent in the case of 
self-adjoint-positive matrices [58]. Moreover, the multilevel algorithm has provably op- 
timal O(N) convergence properties under the standard assumptions for uniform refine- 
ments [94], and is nearly-optimal O(NlogN) under very weak assumptions on adap- 
tively refined problems [12, 2]. In the adaptive setting, a stabilized (approximate wavelet) 
hierarchical basis method is employed [2]. External software can also be used to gener- 
ate the prolongation matrices, so that a number of different graph theory -based algebraic 
multilevel coarsening algorithms may be used to generate the subspace hierarchy. 

Coupled with the superlinear convergence properties of the outer inexact Newton it- 
eration in Algorithm 2.2.2, this leads to an overall complexity of O(N) or 0(iV log N) 
for the solution of the discrete nonlinear problems in Step 1 of Algorithm 2.2.1. Com- 
bining this low-complexity solver with the judicious placement of unknowns only where 
needed due to the error estimation in Step 2 and the subdivision algorithm in Steps 3-6 of 
Algorithm 2.2.1, leads to a very effective low-complexity approximation technique for 
solving a general class of nonlinear elliptic systems on 2- and 3-manifolds. 
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3.5. Parallel computing in MC: The Parallel Partition of Unity Method (PPUM). 

MC incorporates a new approach to the use of parallel computers with adaptive finite 
element methods, based on combining the Partition of Unity Method (PUM) of Babuska 
and Melenk [5] with local error estimate techniques of Xu and Zhou [97]. The algorithm, 
which we refer to as the Parallel Partition of Unity Method (PPUM), is described in detail 
in [9, 13]. The idea of the algorithm is as follows. 

Algorithm: (PPUM - Parallel Partition of Unity Method [9]) 

(1 ) Discretize and solve the problem using a global coarse mesh. 

(2) Compute a posteriori error estimates using the coarse solution, and decompose 
the mesh to achieve equal error using weighted spectral or inertial bisection. 

(3) Give the entire mesh to a collection of processors, where each processor will per- 
form a completely independent solve-estimate-refine loop (Step 2 through Step 6 
in Algorithm 2.2.1), restricting local refinement to only an assigned portion of 
the domain. The portion of the domain assigned to each processor coincides 
with one of the domains produced by spectral bisection with some overlap (pro- 
duced by conformity algorithms, or by explicitly enforcing substantial overlap ). 
When a processor has reached an error tolerance locally, computation stops on 
that processor. 

(4) Combine the independently produced solutions using a partition of unity subor- 
dinate to the overlapping subdomains. 

While the algorithm above seems to ignore the global coupling of the elliptic problem, 
some recent theoretical results [97] support this as provably good, and even optimal in 
some cases. The principle idea underlying the results in [97] is that while elliptic prob- 
lems are globally coupled, this global coupling is essentially a "low-frequency" cou- 
pling, and can be handled on the initial mesh which is much coarser than that required 
for approximation accuracy considerations. This idea has been exploited, for example, 
in [95, 96], and is in fact why the construction of a coarse problem in overlapping do- 
main decomposition methods is the key to obtaining convergence rates which are inde- 
pendent of the number of subdomains (c.f. [94]). A more complete description can be 
found in [9], along with examples using MC and the 2D adaptive finite element package 
PLTMG [10]. An analysis of the global L 2 - and iJ 1 -error in solutions produced by the 
algorithm appears in the next section. An example showing the types of local refinements 
that occur within each subdomain is depicted in Figure 5. 




Figure 5 . An example showing the types of local refinements that are 
created by PPUM. 



3.6. Global L 2 - and i/^-error estimates for PPUM. In order to analyze the error be- 
havior in PPUM, we first review the partition of unity method (PUM) of Babuska and 
Melenk [5]. Let fi C l d be an open set and let be an open cover of with a 
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bounded local overlap property: For all igSI, there exists a constant M such that 

sup{ i | x G Qi } < M. (3.6) 

i 

A Lipschitz partition of unity {(pi} subordinate to the cover {fij} satisfies the following 
five conditions: 

J2<t>i( x ) = 1, Varefi, (3.7) 

% 

(j>i G C k (n) Vi, (A; > 0), (3.8) 

sup 4>i C fij, Vi, (3.9) 

< Coo.Vi, (3.10) 

The partition of unity method (PUM) builds an approximation u ap = ^ </>jt>j where the 
Vi are taken from the local approximation spaces: 

Vi c c k (n n n 4 ) c # x (fi n a), Vi, (fc > o). (3.12) 

The following simple lemma makes possible several useful results. 



Lemma 3.1. Let w, Wi G H l (Q) with supp Wi C £7 n f2j. TTien 
J^IMIn^) < M ||«;||^ (n) , A; = 0, 1 

i 

\\J2 Wi W 2 H k m < MY^\\wi\\ 2 HH(in[k) , k = 0,l 

i i 

Proof The proof follows from (3.6) and (3.7)-(3.1 1); see [5]. □ □ 

The basic approximation properties of PUM are as follows. 

Theorem 3.2 (Babuska and Melenk [5]). If the local spaces Vi have the following ap- 
proximation properties: 

\\u - u»||z2(nnfii) < e Q (i), Vi, 
||V(«-Vi)||i,3(nnno < e i(*)> 
then the following a priori global error estimates hold: 

1/2 

\u — u„ 



1/2 



Proof. This follows from Lemma 3.1 by taking w — u ap = J2i4>i(u — Vi) and u>j = 

4>i(u-Vi). □ □ 

We now give a global iJ 1 -error estimate of the PPUM adaptive algorithm proposed 
in [9]. We can view PPUM as building a PUM approximation u w = ^ 0^ where the 
Vi are taken from the local spaces: 

Vi = XiVf C C k (Vl n ^) C if 1 ^ n a), Vi, (A; > 0), (3.13) 
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where is the characteristic function for f^, and where 

Vf c C k (Q) c H\n), Vi, (Jfe > 0). (3.14) 

In PPUM, the global spaces V^ 9 in (3.13)— (3.14) are built from locally enriching an ini- 
tial coarse global space Vo by locally adapting the finite element mesh on which Vq is 
built. (This is in contrast to classical overlapping schwarz domain decomposition meth- 
ods where local spaces are often built through enrichment of Vq by locally adapting the 
mesh on which Vq is built, and then removing the portions of the mesh exterior to the 
adapted region.) The PUM space V is then 



V = < v I v 



V V 



<fiiVi, Vi e Vi | 

i i J 



Consider now the following linear elliptic problem in the plane: 

-v.(«v«) = /inn 

u = on oil, 

where aij G W^°°(n), f e L 2 {Q), ay&& > a > 0, V& ^ 0, where Q C M 2 is a 
convex polygon. (The results below also hold more generally for classes of two- and 
three-dimensional nonlinear problems.) A weak formulation is: 

Find u G Hl{Q) such that (F(u),v) = 0, Vf G H^Q), 

where 



(F(u),v) = / aV« • Vf dx — / fv dx. 
Jn Jq 

The PUM is usually used to solve a PDE as a Galerkin method in the globally coupled 
PUM space (cf. [46]): 

Find u ap G V C H^Q) s.t. (F(u ap ),v) = 0, Vf G 7 C #o(fi). 

In contrast, PPUM proposed in [9] builds an approximation u pp from decoupled local 
Galerkin solutions: 



^0iWi = ^0iwf, (3.16) 



where each u\ satisfies: 

Find uf G Vf such that (F(v%),vf) = 0, Vff G Vf . (3.17) 

We have the following global error estimate for the approximation u pp in (3.16) built 
from (3.17) using the local PPUM parallel algorithm. 

Theorem 3.3. Assume the solution to (3.15) satisfies u G H 1+a (Q), a > 0, and assume 
that quasi-uniform meshes of sizes h and H > h are used for VL® and VL\VL® respectively. 
If 'diarn(fij) > 1/Q > Vi then the global solution u pp in (3.16) produced by the PPUM 
Algorithm 3.5 satisfies the following global error bounds: 



\u — u 



pp\\L 2 (Q) 



< y/PMCoo (dh a + C 2 H l+a ) , 



V(u~u pp )\\ L 2 {n) < J2PM(Q 2 C 2 G + Cl)(C 1 h a + C 2 H 1+a ), 
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where P = number of local spaces Vi. Further, if H < h a ^ 1+a ^ then: 
\\u - UppWtf^) < V PMC^ maxICx, C 2 }h a , 



\\V(u-u pp )\\ LHn) < ^2PM(QiC 2 G + Cl)m a x{C 1 ,C 2 }h a , 
so that the solution produced by Algorithm 3.5 is of optimal order in the H x -norm. 

Proof. Viewing PPUM as a PUM gives access to the PUM a priori estimates in Theo- 
rem 3.2; these require local estimates of the form: 

\\u - UiHz^nnfii) = \\u - ufHz^nnfii) < e (i), 
||V(w - Ui)\\ L 2 iUnni) = ||V(u - wf)||L2 (nnn .) < e x (i). 

Such local a priori estimates are available for problems of the form (3.15) [74, 97]. They 
can be shown to take the following form: 



\\ u - u i llffH^nn) < C inf \\u - v { \\ H ^n°nn) + \\ u ~ u i\\^(n) 
where 

V° C C fe (fi° nfl)c H\Q, t n fi), 

and where 

Since we assume u E H 1+a (Vl), a > 0, and since quasi-uniform meshes of sizes h and 
H > h are used for and respectively, we have: 

Ik-MfllHH^nn) = ( k ll«-«illLa(n i nn) + ll V ( w - w i) ll£a(n 4 nn)J 
< dh a + C 2 H 1+a . 

I.e., in this setting we can use e (i) = ei(i) = C\h a + C 2 H 1+a . The a priori PUM 
estimates in Theorem 3.2 then become: 

|V(u- 1^)11x2(0) < v/2M 

1/2 



2 



^ ( diam(0^ ) + C ~ 



(d/i 01 + C 2 H 



l+a\2 



If P = number of local spaces and if diam(f2j) > 1/Q > Vi, this is simply: 

1^(0) < VPMC^ (dh a + C 2 H 1+a ) , 



\Mu - u pp )\\ LHn) < \JlPM{QiCl + CI) (C x h a + C 2 H 1+a ) . 

If H < /?, a /( 1+a ) then u pp from PPUM is asymptotically as good as a global Galerkin 
solution when the error is measured in the if 1 -norm. □ □ 

Estimates similar to Theorem 3.3 appear in [97] for a variety of related parallel algo- 
rithms. Note that improving the estimates in the L 2 -norm is not possible; the required 
local estimates simply do not hold. Improving the solution quality in the L 2 -norm would 
require more global information. 
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3.7. Availability of MC and the supporting tools MALOC and SG. MC is built on 
top of a low-level portability library called MALOC (Minimal Abstraction Layer for 
Object-oriented C). Most of the images appearing in this paper were produced using 
a software tool called SG (Socket Graphics), which is also built on top of MALOC. 
MALOC, MC, and SG were developed by the author over several years, with generous 
contributions from a number of colleagues. MALOC, MC, and SG are freely redis- 
tributable under the GNU General Public License (GPL), and the source code for all 
three packages is freely available at the following website: 



MALOC, MC, and SG, as well as a fully functional MATLAB version of MC called 
MCLAB, are part of a larger project called FETK (The Finite Element Toolkit). Infor- 
mation about FETK can be found at: 



4. Example: The Hamiltonian and momentum constraints in the 



The evolution of the gravitational field was conjectured by Einstein to be governed by 
twelve coupled first-order hyperbolic equations for the metric of space-time and its time 
derivative, where the evolution is constrained for all time by a coupled four-component 
elliptic system. This four-component elliptic system consists of a nonlinear scalar Hamil- 
tonian constraint, and a linear 3-vector momentum constraint. The evolution and con- 
straint equations, similar in some respects to Maxwell's equations, are collectively re- 
ferred to as the Einstein equations. Solving the constraint equations numerically, sepa- 
rately or together with the evolution equations, is currently of great interest to the physics 
community (cf. [56, 55, 26] for more detailed discussions of this application). 

The Hamiltonian and momentum constraints in the Einstein equations, taken sepa- 
rately or together as a coupled system, have the form (2.5)-(2.7). Allowing for both 
Dirichlet and Robin boundary conditions as are typically used in black hole and neutron 
star models (cf. [56, 55, 26]), the strong form can be written as: 



http : / / www . scicomp . ucsd . edu/ "mholst / 



http : / / www . f etk . org 



Einstein equations 



A0 



1 



£0+l(trK) 2 5 

l(*A ab + {LW) ab ) 2 <p- 7 - 27TP0- 3 in M 



(4.1) 



8 



n a D a <j) + c<j) = zondiM, 
4> = fondoM, 



(4.2) 
(4.3) 



D b {LW) ah = -4) & b a tiK + 8vrf in M 
3 

(LW) ab h b + C\W b = Z a ond 1 M, 
W a = F a ond M, 

where the following standard notation has been employed: 



(4.4) 



(4.5) 
(4.6) 



A<f) = D a D a <f), 
(LW) ab = D a W b + b b w a - -i ab D c W c , 



trK = 7 ab K ab , 

I ft \2 s-iabs-i 
[yah) — t> ^ab- 
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The symbols in the equations (R, K, *A ab , p, j a , z, Z a , f, F a , c, and Cg) represent various 
physical parameters, and are described in detail in [56, 55, 26] and the referenences 
therein. 

Equations (4.1)-(4.6) are known to be well-posed only for restricted problem data 
and manifold topologies [77, 75, 76]. Below we will present two well-posedness results 
from [56] which hold under certain assumptions. Note that if multiple solutions in the 
form of folds or bifurcations are present in solutions of (4.1)-(4.6) then path-following 
numerical methods will be required for numerical solution [62, 63]. 

4. 1 . Weak formulation, linearization, and well-posedness. Both the Hamiltonian con- 
straint (4.1) and the momentum constraint (4.4), taken separately or as a system, fall into 
the class of second-order divergence-form elliptic systems of tensor equations in (2.5)- 
(2.7). Derivation of the weak formulation produces a weak system of the form (2.14)- 
(2.15), with some interesting twists along the way described in [56]. Following the no- 
tation in [56], we employ a background (or conformal) metric to define the volume 
element dx = \/det % b dx 1 dx 2 dx 3 and the corresponding boundary volume element ds, 
and for use as the manifold connection for covariant differentiation. The notation for 
covariant differentiation using the conformal connection will be denoted D a to be con- 
sistent with the relativity literature, and the various quantities from Section 2.1 will now 
be hatted to denote use of this conformal metric. For example, the unit normal to dAA 
will now be denoted h a . 

Ordering the Hamiltonian constraint first in the system (2.5), and defining the product 
metric Qij and the vectors u l and v J appearing in (2.1 1) and (2.15) as: 



1 

g a b 



u 



W a 



v b 



it is shown in [56] that the coupled Hamiltonian and momentum constraints have a cou 
pled weak formulation in the form of (2.14), where the form definition is as follows: 

(F(u),v) = (F([<p,W a ]),[i/;,V a }) = (F H (4>),tP) + (F M (W a ),V a ). 

The individual Hamiltonian form is shown in [56] to be: 



(4.7) 



D a (j)D a ^ dx + / P'((j))ipdx 



M 



(4.8) 



M 



where 



+ / (c(j) — z)ij) ds, 

'diM 



P'{cf>) = X -R4> + ^(tiK) 2 ^ - ^(*A ab + (LW) ab ) 2 0- 7 - 2np4>- 3 , 



and the momentum form is shown in [56] to be: 



(F M (W a ),V a ) 
'2 



M 



2p(EW) ab (EV) ab + XD a W a D b V b ) dx 



(4.9) 



(4.10) 



6 D a trK + 87rj a ) V a dx I {C a b W b - Z a )V a ds, 

where p = 1, X = —2/3, and where the deformation tensor (EV) ab is the symmetrized 
gradient: 

(EV) ab = - (t) b V a + D a V b ) . (4.11) 
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The Gateaux-derivative of the nonlinear weak form (F(-), ■) in equation (4.7) above 
is needed for use in Newton-like iterative solution methods such as Algorithm 2.2.2. 
Defining an arbitrary variation direction «; = [£, X a ], it is shown in [56] that the Gateaux- 
derivative takes the following form (for fixed [0, W a }), linear separately in each of the 
variables [f , X a ] and [if>, V a \: 

(DF{[<l>,W*\)[tL,X%ty),V a ])= (4.12) 
/ (c^ + C%X b V a ) ds 

+ J (p a £b a il) + 2p(EX) ab (EV) ab + \b a x a b b v h } dx 
+ J ^( tri W + \(*A ab + {LWU) 2 ^ 8 + foMrA & dx 

- (^{*A ab + {LW) ab )r 7 ^j (LX) ab tP dx + (^L^trKj V a i dx. 

Now that the nonlinear weak form (F(-), •) and the associated bilinear linearization form 
(DF(-)-, ■) are defined and can be evaluated using numerical quadrature, the assembly of 
the nonlinear residual as well as linearizations about any point can be performed precisely 
as outlined above for a generic nonlinear finite element method. Again, once we have the 
weak formulation and a linearization, the discretization in MC is automatic and generic. 
Although the forms (F(-), ■) and (DF(-)-, ■) above appear somewhat complicated in the 
case of the constraint equations, the discretization in MC involves simply evaluating the 
integrands for use in quadrature formulae. 

We now state two new existence and uniqueness results for the Hamiltonian and mo- 
mentum constraints which were established recently in [56]. A number of assumptions 
on the problem data are required. 

Assumption 4.1. We assume that Ai is a connected compact Riemannian 3-manifold 
with Lip schitz- continuous boundary dAA. We also assume that the data has the following 
properties: 

lab e W X >°°{M), K ab e W lfil \M), <p e L°°(M), W a e H\M), 

f e H-\M), C\ e L 2 (9i7W), Z a e L^{d x M), F a e H 1/2 (d M), 
where W \d Q M — F a in the trace sense, and where for some constant a > 0, 

[ C\V h V a dx > a\\V a \\ 2 LHdiM) , V V a E L 4 (diM). 

Jd x M 

Assumption 4.2. We assume that Ai is a connected compact Riemannian 3-manifold 
with Lipschitz-continuous boundary dAi, where meas(<9o-M) > 0. We also assume that 
the data has the following properties: 

i ab e W l >°°(M), K ab e L°°{M), W a e W^°°{M), 4> e H\M) n L°°(M), 

R,*A ab , P eL°°(M), czeL^id.M), f e # 1/2 (d -M) n L°°(d M), 
where 4>\d M — f in the trace sense, and 

< inf / < / < sup / < oo, a.e. in doAi, 

d oM q qM 

p > 0, a.e. in Ai, c > 0, z = 0, a.e. on d\AA. 
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Theorem 4.3. Let Assumption 4.1 hold. Then there exists a unique solution W a G W a + 
Hq D (M) to the momentum constraint equation (4.4)-(4.6) which depends continuously 
on the problem data. Moreover, U a = W a — W a G D (M) satisfies the following a 
priori bound: 

\\U a \\m(M) < ||^ a |U 2 (>!) + — , (4.13) 

where a is the strong ellipticity constant and L is a bound on the linear functional arising 
in the weak form. 7/"meas(9 A^) > 0, then the following bound also holds: 

\\U a \\m { M) < -, (4-14) 
m 

where m is the coercivity constant. 

Proof. The proof given in [56] is based on the use of a Riesz-Schauder alternative ar- 
gument (uniqueness implies existence), which is accessible after establishing that the 
momentum weak form operator has a number of properties, including strong ellipticity 
and satisfaction of a Garding inequality. □ □ 

Theorem 4.4. Let Assumption 4.2 hold. Then there exists a unique solution <p G + 
Hq D (M) to the Hamiltonian constraint equation (4.1)— (4.3). The solution satisfies a 
priori L°° -bounds and is strictly positive a.e. in M.. 

Proof. The proof given in [56] is based on variational analysis and fixed-point argu- 
ments, after using a weak maximum principle to remove the poles at the origin in the 
nonlinearity. □ □ 

The two results above indicate that the momentum and Hamiltonian constraints on 
connected compact Riemannian manifolds with Lipschitz boundaries have well-posed 
weak formulations in the unweighted Sobolev spaces Hq D (JA). A small amount of 
additional regularity, namely the intersection of Assumptions 4.1 and 4.2, is required 
to give simultaneously well-posed weak formulations. Under smoothness assumptions 
on the boundary and coefficients, this minimal additional regularity can be shown for 
both <p and W a using elliptic regularity arguments (cf. [70] for a discussion of the linear 
elasticity case which can be adapted here for the momentum constraint). Unfortunately, 
elliptic systems such as the momentum constraint do not satisfy maximum principles 
analogous to the weak maximum principle derived for the (scalar) Hamiltonian constraint 
in [56], and as a result it is more difficult to establish L°° -bounds on W a . Note that 
simultaneous well-posedness of the Hamiltonian and momentum constraints individually 
does not imply well-posedness of the coupled system. Limited results for the coupled 
system exist for some simplified situations; cf. [60, 59, 34, 33, 77]. Some new results for 
the coupled system, based on Theorems 4.3 and 4.4, appear in [56]. 

4.2. Quasi-optimal a priori error estimates for Galerkin approximations. In this 
section we consider the theory for Galerkin approximations of the Hamiltonian and mo- 
mentum constraints. Following the approaches in [29, 84, 61] for related problems, we 
establish two abstract results, the first of which applies to linear variational problems sat- 
isfying a Garding inequality, whereas the second result applies to monotonically nonlin- 
ear variational problems. When applied to the Hamiltonian and momentum constraints, 
each result will take the form: 



(4.15) 
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where Uh is a Galerkin approximation such as provided by a finite element discretiza- 
tion, and where Vh C Hq D (A4) is the subspace of continuous piecewise polynomials 
defined over simplices. These results are quasi-optimal in the sense that they imply that 
a Galerkin solution of either the Hamiltonian or momentum constraint is within a con- 
stant of being the best approximation in the particular subspace in which the Galerkin 
solution lives. After giving the two abstract results along with their simple short proofs, 
we indicate how they can be applied to the momentum and Hamiltonian constraints in 
the context of Galerkin finite element methods. 

While the term on the left in (4.15) is in general difficult to analyze, the term on the 
right represents the fundamental question addressed by classical approximation theory in 
normed spaces, of which much is known. To bound the term on the right from above, one 
picks a function in Vh which is particularly easy to work with, namely a nodal or general- 
ized interpolant of it, and then one employs standard techniques in interpolation theory. 
Therefore, it is clear that the importance of approximation results such as (4.15) are that 
they completely separate the details of the momentum and Hamiltonian constraints from 
the approximation theory, making available all known results on finite element interpo- 
lation of functions in Sobolev spaces (cf. [35]). There are some additional difficulties 
in using the standard finite element interpolation theory associated with the fact that we 
are working with a domain with the structure of a Riemannian 3 -manifold rather than 
an open set in M d ; these are being addressed in work in progress [53], and will not be 
discussed in detail here. 

4.2. 1 . Approximation theory for the momentum constraint. We now give a quasi-optimal 
a priori error estimate which characterizes the quality of a Galerkin approximation to 
the solution of the momentum constraint. Quasi-optimal estimates are quite standard 
in the finite element approximation theory literature for V-elliptic bilinear forms, but 
unfortunately it is shown in [56] that the momentum constraint weak form is only V- 
coercive (satisfying a Garding inequality). However, following Schatz [84] we show 
this is sufficient to establish similar quasi-optimal results for the momentum constraint 
(cf. [85, 97] for related results). 

In order to derive such a result following the approach in [84], we begin with a Gelfand 
triple of Hilbert spaces V C H = H* C V* with continuous embedding, meaning that 
the pivot space H and its dual space H* are identified through the Riesz representation 
theorem, and that the embedding V C H is continuous. A consequence of this is: 

\\u\\ H < C\\u\\ v , Vit G V, (4.16) 

where we will assume that the embedding constant C = 1 (the norm || • \\y can be rede- 
fined as necessary). In our setting of the momentum constraint, we have H = L 2 (A4) 
and V = Hq d {M) generating the triple; we will stay with the abstract notation involving 
H and V for clarity. We are given the following variational problem: 

Find w G V s.t. A{u,v) = F(v), G V, (4.17) 

where the bilinear form A(u, v) : V x V h->- R is bounded 

A{u,v) < M\\u\\ v \\v\\ v , \/u,v G V, (4.18) 

and V-coercive (satisfying a Garding inequality): 

m||w||y < + A(u, u), Vm G V, where m > 0, (4.19) 

and where the linear functional F(v) : V i— > M. is bounded and thus lies in the dual space 

y*. 

F(v) < L\\v\\ v , G V. 
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It is shown in [56] that the weak formulation of the momentum constraint (4.10) fits into 
this framework; to simplify the discussion, we have assumed that any Dirichlet function 
u has been absorbed into the linear functional F(v) in the obvious way. Our discussion 
can be easily modified to include approximation of u by Uh- 

Now, we are interested in the quality of a Galerkin approximation: 

Find u h G V h C V s.t. A(u h , v) = A(u, v) = F(v), \/v G V h C V. (4.20) 

We will assume that there exists a sequence of approximation subspaces Vh C V param- 
eterized by h, with V/^ C 14 2 when h 2 < hi, and that there exists a sequence {ah}, with 
lim^o a h = 0, such that 

II'" — Uh\\n < CLh\\u — Uh\\v, when A(u — Uh, v ) = 0, Wv G Vh C V. (4.21) 

The assumption (4.21) is very natural; in our setting, it is the assumption that the error in 
the approximation converges to zero more quickly in the L 2 -norm than in the i^-norm. 
This is easily verified in the setting of piecewise polynomial approximation spaces, under 
very mild smoothness requirements on the solution u; cf. Lemmas 2.1 and 2.2 in [97]. 
Under these assumptions, we have the following a priori error estimate. Although the 
assumptions are slightly different, the result and the main idea for the simple proof we 
give below (included for completeness) go back to Schatz [84] (see also [85, 97]). 

Theorem 4.5. Let V C H c V* be a Gelfand triple ofHilbert spaces with continuous 
embedding. Assume that (4.17) is uniquely solvable, and that assumptions (4.16), (4.18), 
(4.19), and (4.21 ) hold. Then for h sufficiently small, there exists a unique approximation 
Uh satisfying (4.20), for which the following quasi-optimal a priori error bounds hold: 

\W — Uh\\v < C hif \\u — v\\v, (4.22) 

vev h 

\\u — Uh\\H < Ca>h inf \\u — v\\v, (4.23) 

vev h 

where C is a constant independent ofh. If K < in (4.19), then the above holds for all 
h. 

Proof. The following proof follows the idea in [84]. We begin with the Garding inequal- 
ity (4.19) and then employ (4.18): 

m\\u - u h \\y - K\\u - u h \\ 2 H < A(u-u h ,u-u h ) 

= A(u — Uh,u — v) 

< M||u-u ft ||v||u-i;||v, (4.24) 

where we have used Galerkin orthogonality: A(u — u h ,v) = 0, Vt> G Vh, to replace Uh 
with an arbitrary v G V h . Excluding first the case that \\u — u h \\ v = we divide through 
by m\\u — u h \\ v and employ (4.16) and (4.21), giving Vt> G V h , 

Kcih\ | K\\u — Uh\\ \ 

\U - Uh\\v < \\U - Uh \\v - 



my m\\u — Uh\\v 

M„ 

< — \\u-v \\ v , (4.25) 
m 

which we note also holds when \\u — Uh\\ =0. 

Assume first that K > 0. Since lim^^o a h = 0, there exists h such that ah < 
mj K, V7i < h. This implies Vt> G Vh, 

( %\ ,, ( Ka h \ n . M„ 

1 * u - u h v < 1 u - u h v < — \ u - v v . (4.26) 

V m I \ m I m 
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Taking u = in (4.20) together with v = in (4.26), with h < h, implies that the 
homogeneous problem 

Find u h G V h s.t. A(u h , v ) = 0, VwG V h , 

has only the trivial solution, so that by the discrete Fredholm alternative a solution Uh 
to (4.20) is unique and therefore exists. Equation (4.26) then finally gives (4.22) when- 
ever h < h, with the choice 

M M 



Assume now that K < 0. Directly from (4.25) we can conclude (4.22) with C = 
M/m, which is completely independent of n; this then becomes Cea's Lemma for V- 
elliptic forms [35]. Moreover, the continuous and discrete problems are both uniquely 
solvable due to V-ellipticity (4.19), independent of h. 

In either case of K > or K < 0, the second estimate (4.23) now follows immediately 
from assumption (4.21). □ □ 

In the case of the momentum constraint it was established in [56] that the assump- 
tions required for Theorem 4.5 hold, with the exception of (4.21). In the case of Robin 
boundary conditions, it was shown in [56] that 1 < a = K < 4/3. This gives 

c= M < M 



«(! ~ a h) 1 - a h 

Under the mild assumption that the a priori bound (4.13) or (4.14) can be shown to hold 
in a slightly stronger Sobolev norm, referred to as an elliptic regularity estimate: 

\\W a \\m+s {M) < — , s > 0, 
m 

then it can be shown that (4.21) holds in the setting of piecewise linear finite element 
spaces, with = rC" 1 for some 7 > 0, where n = dim(Vfo). This makes it clear that 
the requirement that h be sufficiently small is not a practical restriction on applying the 
finite element method to the momentum constraint. 

4.2.2. Approximation theory for the Hamiltonian constraint. We consider now the non- 
linear Hamiltonian constraint, and derive a quasi-optimal a priori error estimate for 
Galerkin approximations analogous to that derived in the previous section for the mo- 
mentum constraint. The approximation theory for Galerkin approximations to the non- 
linear Hamiltonian constraint (4.8) is somewhat more complex than for the momentum 
constraint (4.10). However, it is still possible to establish a result for the Hamiltonian 
constraint which shows that a Galerkin approximation is quasi-optimal under some weak 
assumptions on the nonlinearity. A number of such estimates have appeared in the liter- 
ature; the result we derive below is similar to estimates in [35, 29, 61]. 

We begin again with a Gelfand triple of Hilbert spaces V C H = H* C V* with 
continuous embedding, so that again (4.16) holds. We are given the following nonlinear 
variational problem: 

Find u e V s.t. A(u, v) + (B(u),v) = F(v), \/v E V, (4.27) 

where the bilinear form A(u, v ) : V x V 1— > M. is bounded 

A(u,v) < M\\u\\ v \\v\\v, Vu,v G V, (4.28) 

and V-elliptic: 

571 1 Ml v- — A{u,u), Vw G V, where m > 0, (4.29) 
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where the linear functional F[v) : V (-)■ K is bounded and thus lies in the dual space V*: 

F(v) < L\\v\\ v , G V, 

and where the nonlinear form (B(u),v) : V x V h->- R is assumed to be monotonic: 

< (B(u) — B(v), u — v), Vw, ti £ V, (4.30) 

where we have used the notation: 

(B(u) - B(v),w) = (B(u),w) - (B(v),w). (4.31) 

We are interested in the quality of a Galerkin approximation: 

Find u h G V h s.t. A(u h ,v) + (B(u h ),v) = F(v), G VJ,, (4.32) 

where C V. We will assume that (B(u),v) is bounded in the following weak sense: 
If u G V satisfies (4.27), if Uh G Vh satisfies (4.32), and if v G Vh, then there exists a 
constant K > such that: 

(B(u) — B{uh)i u — v) < K\\u — Uh\\v\\u — v\\v- (4.33) 

It is shown in [56] that the weak formulation of the Hamiltonian constraint (4.8) fits 
precisely into this framework with the possible exception of (4.33); we will show below 
that a priori bounds such as those established in [56] can be used to establish (4.33). 
We have again assumed that any Dirichlet function u has been absorbed into the various 
forms in the obvious way to simplify the discussion. The discussion can be modified to 
include approximation of uby Uh- 

Again, we are interested in the quality of a Galerkin approximation u h satisfying (4.32), 
or equivalently: 

A(u - u h ,v) + (B(u) - B(u h ),v) = 0,\/veV h C V. 

As before, we will assume that there exists a sequence of approximation subspaces Vh C 
V parameterized by h, with C Vh 2 when h 2 < hi, and that there exists a sequence 
{a h }, with lim^o a h = 0, such that 

\\u-u h \\ H <a h \\u-u h \\ v , (4.34) 

holds whenever u h satisfies (4.32). The assumption (4.34) is again very natural; see the 
discussion above following (4.21). Under these assumptions, we have the following a 
priori error estimate. 

Theorem 4.6. Let V C H C V* be a Gelfand triple of Hilbert spaces with contin- 
uous embedding. Assume that (4.27) and (4.32) are uniquely solvable, and that as- 
sumptions (4.16), (4.28), (4.29), (4.33), and (4.34) hold. Then the approximation Uh 
satisfying (4.32) obeys the following quasi-optimal a priori error bounds: 

\\ u — u h\\v < C inf \\u — v\\v, (4.35) 

vev h 

\\ u — uii\\h < Ca>h inf Ilit — v\\v, (4.36) 

vev h 

where C is a constant independent of h. 

Proof. We begin by subtracting (4.32) from (4.27), and taking v — w E Vh CV in both 
equations, giving: 

A(u - u h , w) + (B(u) - B(u h ),w) = 0, Ww G V h . (4.37) 
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In particular if v E Vh, so that w = v — u h E V h , this implies that 

A(u - u h ,v - u h ) = (B(u h ) - B(u),v - u h ) 

= (B(u h ) - B(u),v - u) - (B{u h ) - B(u),u h - u) 

< {B{u h )-B{u),v-u), (4.38) 

where we have employed monotonicity (4.30). Beginning now with (4.29) we have for 
arbitrary v E V h that 

m\\u-u h \\y < A(u-u h ,u-Uh) 

= A(u-u h ,u-v) + A(u-u h ,v -u h ) 

< A(u — Uh,u — v) + (B(uh) — B(u),v — u) 

< M||u-tt fc ||v||u-u||v (4.39) 

— lt/j||y||u — v\\y. 

where we have used (4.38), (4.28), and (4.33). Excluding first the case that ||m — w?i ||y — 
we divide through by m\\u — Uh\\v, giving 

/ M + K\ „ 

u - u h \\ v < u - v\\ v , G V h , (4.40) 

\ m J 

which we note also holds when \\u — Uh\\ = 0. This gives (4.35) with C = (M + K)/m. 
The second estimate (4.36) now follows immediately from assumption (4.34). □ □ 

In the case of the Hamiltonian constraint the nonlinear weak form (B(u),v) has the 
form 

(B{u),v) = I P'{u)vdx, 
Jm 

where P'(u) is defined in (4.9). If both u and Uh satisfy a priori bounds as established 
in [56], then the continuity of P"(x) on (0, oo) implies that there exists w E L°°(A4) 
satisfying similar bounds such that 

P\u) - P'(u h ) = P"(w)(u - u h ), a.e. in M. 

Consider now 

(B{u)-B{u h ),u-v) 
{P'{u) - P'(u h )){u-v) dx 

M 

P"(w)(u — Uh){u — v) dx 

M 

< \\P"{w)\\l™{<x<w<I3)\\u - U h \\ L 2 {M) \\u - v\\ L 2( M) 

< K\\u - u h \\ H i {M) \\u - v\\ H i {M) . 

Therefore, (4.33) holds with K = \\P"(w) \\L°°(a<w<p), which can be computed explicitly 
from the results in [56]. Although a Galerkin approximation Uh constructed from finite 
element bases will not in general satisfy a discrete maximum principle which would lead 
to a priori bounds as in [56], it is possible to establish L°° -bounds for a Galerkin finite 
element solution to the Hamiltonian constraint under some assumptions on the size and 
shape of the elements in the mesh (cf. Theorem 3.2 in [64]). 

Therefore, we see that in the case of the Hamiltonian constraint we have established 
that the assumptions required for Theorem 4.6 hold, with the exception of (4.34). Under 
the mild additional regularity assumption: 

\m+s(M) <C<oo, s > 0, 
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where C depends on the data, then it can be shown that (4.34) holds in the setting of 
piecewise linear finite element spaces, with = n~ 7 for some 7 > 0, where n = 
dim(V h ). 

Other approaches also lead to well-posed weak formulations of the Hamiltonian con- 
straint with associated approximation theory. In particular, an obstacle problem for- 
mulation is possible as a technique for handling the pole at the origin in the Hamil- 
tonian constraint, leading to a nonlinear variational inequality. This approach requires 
fewer assumptions on the data in the Hamiltonian constraint than we have assumed 
here. Although several difficulties arise in a priori error analysis, a number of results 
for linear and nonlinear variational inequalities are known, and could be applied in this 
case. Results similar to Theorem 4.6 are obtainable under the same minimal assump- 
tion G H l (M) required to give a well-posed weak formulation (cf. [35, 29, 30, 44]). 
Introducing a cut-off function in place of the two terms with poles in the Hamiltonian 
constraint leads to a well-posed weak formulation, although the error analysis is not 
clear. Approaches based on weighted Sobolev spaces also lead to well-posed weak for- 
mulations, but incorporation of weights into the finite element subspaces is technically 
complicated. 

4.3. Numerical solution using MC. To use MC to calculate the initial bending of space 
and time around two massive black holes separated by a fixed distance by solving the 
above constraint equations, we place two spherical objects in space, the first object hav- 
ing unit radius (after appropriate normalization), the second object having radius 2, sep- 
arated by a distance of 20. Infinite space is truncated with an enclosing sphere of radius 
100. (This outer boundary may be moved further from the objects to improve the accu- 
racy of boundary condition approximations.) Resonable choices for the remaining func- 
tions and parameters appearing in the equations are used below to completely specify the 
problem for use as an illustrative numerical example. (More careful examination of the 
various functions and parameters appear in [56], and a number of detailed experiments 
with more physically meaningful data appear in [55, 26].) 

We then generate an initial (coarse) mesh of tetrahedra inside the enclosing sphere, 
exterior to the two spherical objects within the enclosing sphere. The mesh is generated 
by adaptively bisecting an initial mesh consisting of an icosahedron volume filled with 
tetrahedra. The bisection procedure simply bisects any tetrahedron which touches the 
surface of one of the small spherical objects. When a reasonable approximation to the 
surface of the spheres is obtained, the tetrahedra completely inside the small spherical 
objects are removed, and the points forming the surfaces of the small spherical objects 
are projected to the spherical surfaces exactly. This projection involves solving a lin- 
ear elasticity problem, together with the use of a shape-optimization-based smoothing 
procedure. The smoothing procedure locally optimizes the shape measure function in 
equation (3.1) for a given d-simplex s, in an iterative fashion. A much improved binary 
black hole mesh generator has been developed by D. Bernstein; the new mesh generator 
is described in [55, 26] along with a number of more detailed examples using MC. 

The initial coarse mesh in Figures 6-8, generated using the procedure described above, 
has approximately 31,000 tetrahedral elements and 6,000 vertices. To solve the problem 
on a 4-processor computing cluster using PPUM (see Section 3.5), we begin by partition- 
ing the domain into four subdomains (shown in Figures 9-10) with approximately equal 
error using the recursive spectral bisection algorithm described in [9]. The four subdo- 
main problems are then solved independently by MC, starting from the complete coarse 
mesh and coarse mesh solution. The mesh is adaptively refined in each subdomain until 
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a mesh with roughly 50000 vertices is obtained (yielding subdomains with about 250000 
simplices each). 

The resulting refined subdomain meshes are shown in Figures 11-12. The refinement 
performed by MC is confined primarily to the given region as driven by the weighted 
residual error indicator from Section 2.3, with some refinement into adjacent regions 
due to the closure algorithm which maintains conformity and shape regularity. The four 
problems are solved completely independently by the sequential adaptive software pack- 
age MC. One component of the solution (the conformal factor 0) of the elliptic system 
is depicted in Figure 13 (the subdomain zero solution) and in Figure 14 (the subdomain 
two solution). 

While this example illustrates some of the capabilities of MC, a number of more de- 
tailed examples involving the contraints, using more physically meaningful data, appear 
in [55, 26]. 




Figure 6. The coarse binary black hole mesh (approximately 6,000 ver- 
tices and 31,000 simplices). 



5. Summary 

In this paper we considered the design of adaptive multilevel finite element methods 
for certain elliptic systems arising in geometric analysis and general relativity. We be- 
gan with a brief introduction to nonlinear elliptic tensor systems on manifolds, and then 
discussed adaptive finite element methods for this class of problems. We derived two a 
posteriori error indicators, one of which was local residual-based, and one of which was 
based on a global linearized adjoint or dual problem. 

The implementation of these methods and indicators in the ANSI C finite element 
software package MC was discussed, including detailed descriptions of some of the 
more interesting algorithms and data structures it employs. MC was designed by the 
author specifically for solving general second-order nonlinear elliptic systems of tensor 
equations on Riemannian manifolds with boundary, including domains requiring multi- 
ple coordinate systems. The key feature of MC which makes it particularly useful for 
highly complex tensor systems of PDEs arising in geometric analysis and general rela- 
tivity is its abstraction; in addition to the support for multi-chart manifolds, the MC user 
supplies only two ANSI C functions representing the weak form of the tensor system 
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Figure 7 . Exploded view of the coarse binary black hole mesh showing 
the two interior hole boundaries. 
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Figure 8. Closeup of the interior of the coarse binary black hole mesh. 
The interior holes surfaces of black hole coarse mesh; the larger hole 
surface is colored yellow, the smaller hole surface is colored purple, and 
the exterior boundary is colored red. 

(F(u),v) along with its linearization form (DF(u)w,v). Moreover, the forms them- 
selves may be implemented almost exactly as they are written on paper, due to the fact 
that the quadrature-based assembly allows for tensor expressions to be treated discretely 
as point tensors rather than tensor fields. If residual-based or duality -based a posteriori 
error estimation is to be used, then the user must provide a third function F(u), which 
is essentially the strong form of the differential equation as needed for the residual and 
duality indicators given in Section 2.3. We also described an unusual approach taken 
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Figure 9. Subdomains 2 (red) and 4 (yellow) from spectral bisection of 
the coarse binary black hole mesh; these subdomains enclose two smaller 
subdomains that contain the inner holes. 




Figure 10. Subdomains 3 (blue) and 1 (green) from spectral bisection 
of the coarse binary black hole mesh; these subdomains each contain one 
of the inner holes. 




Figure 1 1 . Closeup of the subdomain 1 refined mesh around the sur- 
face of the smaller hole. (Approximately 51,000 vertices and 266,000 
simplices; only faces of tetrahedra on the boundary surfaces are shown). 
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Figure 12. Closeup of the subdomain 3 refined mesh around the surface 
of the larger hole. (Approximately 45,000 vertices and 228,000 simplices; 
only faces of tetrahedra on the boundary surfaces are show). 




Figure 1 3 . The conformal factor <p from the adapted subdomain 1 solve. 

in MC for using parallel computers in an adaptive setting, based on joint work with R. 
Bank [9]. We then derived global L 2 - and iJ 1 -error estimates for the solutions produced 
by the parallel algorithm, by interpreting the algorithm as a special partition of unity 
method [5] and by using the recent local estimates of Xu and Zhou [97]. 

As an illustrative example, we took a brief look at the Hamiltonian and momentum 
constraints in the Einstein equations. We first summarized a number of operator prop- 
erties and solvability results recently established in [56], and then derived two a priori 
error estimates for Galerkin approximations, completing the theoretical framework for 
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Figure 14. The conformal factor <p from the adapted subdomain 3 solve. 

effective use of adaptive multilevel finite element methods. We finished by presenting 
an illustrative example using the MC software. More detailed examples may be found 
in [55, 26]. 
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